•.*.  RtPOK*  StCjRlif  wokSSif'CATiON 

»cc-' 


;*.  S£C-R.TY  C^SSl^lCATlO^  aht^ob.tv 


AD-A209  018 


64.  ^AM;  C‘  PtRfORMING  organization 

Massachusetts  Institute  of 

Technology,  Civil  Engineerln 


6c  aOORESS  (Oty.  Swfe,  snd  ZIP  Coat} 

77  Massachusetts  Avenue  ,  Room  1-175 
Cambridge,  MA  02129 


REPORT  DOCUMENTATION  PAGE 


10.  RtSTRiCTiVl  MARAINGS 


3.  DiSTRIBJTiON,  awa.,-ABiJTY  Of  REPORT 

Approved  for  public  release; 
distribution  unlimited. 


5.  MONITORING  organization  REPORT  NjMBERiS) 

OiS^iJLO  .  6^-01  ^ 


7a.  NAME  Of  MONITORING  ORGANIZATION 

U.  S.  Army  Research  Office 


7b.  ADDRESS  (Oty,  Sure,  and  ZIP  Coot) 

P.  0.  Box  12211 

Research  Triangle  Park,  NC  27709-2211 


bi.  name  of  funding /sponsoring 
organization 

V.  S.  Arr.v  Research  Office 


8b.  office  symbol  9.  PROCUREMENT  INSTRUMENT  IDENTIFICATION  NUMBER 
(If  spplicsbit) 

sDf)fiL0  2'^'?''/Z-0  0  0  S' 


ADDRESS  (City.  Sufe,  ana  ZIP  Coot) 

10  SOURCE  0=  FUNDING  NUMBERS 

\  0.  Lex  1:211 

.esearch  Triangle  Park,  NC  27709-2211 

PROGRAM 
ELEMENT  NO. 

PROJE'CT 

NO. 

TAS< 

NO. 

T  (incluae  Secur::y  CussmcttionJ 

Mechanics  of  Damage  in  Rate-Sensitive  Construction  Materials 


12  PERSONAL  AUTHOR(S) 


Shyam  Sunder,  V.;  Wu,  Maos 


■ii.  type  of  report 
Technical 


••£.  SUPPLEMENTARY  NOTATION 


j'.3b.  TIME  COVERED 

K.  Date  O'  REPORT  (Year.  Month.  Dty) 

■!  II  iM 

December  30,  1987 

The  viev,  opinions  and/or  findings  contained  in  this  report  are  those 
author (s)  .and  should  not  .be .  construed  as,  an  official  Deoartment  of  the  Arrv  position 


COSATI  CODES 


GROUP  I  SUB-GROUP 


18.  SUBJECT  TERMS  (Continut  on  rtvtnt  if  ntctaary  and  identify  by  block  numotr) 
Damage  mechanics;  rate-sensitive  construction  materials; 
Damage  models;  advanced  construction  materials 


"S  aBSTPaCT  (Continue  on  rtvtrtt  if  ntctssary  and  latnxify  by  block  numotr) 

The  objective  of  this  research  is  to  develop  an  understanding  of  and  characterize 
Che  rate  and  damage  processes  and  their  interaction  in  advanced  construction  materials. 
This  report  presents:  a  multiaxial  differential  flow  model;  a  rate  and  pressure  ^ 
sensitive  anisotropic  damage  model;  and,  specification  for  testing  materials.  j  j 


OTIC" 

f-  LECTE  0^ 

jUN  1  6  1989 


I 


::  D'STRisuTiON/ availability  qp  abstract 
□  UNCAGSIFIED/UNLIMITED  □  SAME  AS  ROT.  □  DTIC  USERS 

21.  abstract  SECURITY  C-ASSIFICaTION 

Unclassified 

NAME  OF  RESPONSIBLE  INDIVIDUAL 

ZZo.  TELEPHONE  (/ncJuoc  Area  CoOt) 

ZZc.  OFFICE  symbol 

CD  FORM*. <73, 8A MAR 


83  APR  eolXion  m«y  d«  used  ur.tii  ex*.*ur;ed. 
All  otner  editions  arc  oosoicie. 


SECURITY  classification  Qt  TH'S  FACE 

UNCL/.SSiriZD 


89  6  16  025 


MECHANICS  OF  DAMAGE 


IN  RATE-SENSITIVE  CONSTRUCTION  MATERIALS 


TECHNICAL  REPORT 


By 


Professor  S.  Shyam  Sunder 
Associate  Professor  of  Civil  Engineering 


Mao  S.  Wu 

Graduate  Research  Assistant 


DECEMBER  30,  1937 

U.S.  ARMY  RESEARCH  OFFICE 
Contract  No.  DAAL03-87-K-0005 


Aceessloa  For 

NTIS  GRAtl 

DTIC  TAB 

Unannounced 
Justification _ 

□ 

By 

Distribution/ 


Availability  Codes 
Avai'l  and/or 
Special 


MASSACHUSETTS  INSTITUTE  OF  TECHNOLOGY 


APPROVED  FOR  PUBLIC  RELEASE; 
DISTRIBUTION  UNLIMITED 


THE  VIEWS,  OPINIONS,  AND/OH  FINDINGS  CONTAINED  IN  THIS  REPORT  ARE 
THOSE  OF  THE  AUTHORS  AND  SHOULD  NOT  BE  CONSTRUED  AS  AN  OFFICIAL 
DEPARTMENT  OF  THE  ARMY  POSITION,  POLICY,  OR  DECISION,  UNLESS  SO 
DESIGNATED  BY  OTHER  DOCUMENTATION. 


I  • 

1 


1 .  TABLE  OF  CONTENTS 


Page 

1.  TABLE  OF  CONTENTS .  1 

2.  LIST  OF  APPENDICES .  2 

3.  TECHNICAL  REPORT .  3 

A.  Statement  cf  Problem .  3 

Introduction .  3 

Objectives  of  Research .  4 

B.  Summary  of  Research  Accomplishments .  4 

Multiaxial  Differential  Flow  Model .  5 

Rate-Sensitive  Continuum  Damage  Model .  11 

Experimental  Facilities  and  Program .  17 

C.  List  of  Publications .  19 

D.  List  of  Participating  Scientific  Personnel .  19 

4.  BIBLIOGRAPHY .  20 


5.  APPENDICES 

Appendix  A 
Appendix  B 


22 


I 


Appendix  A 
Appendix  B 


2 


2.  LIST  OF  APPENDICES 


-  A  MULTIAXIAL  DIFFERENTIAL  FLOW  LAW  FOR 
POLYCRYSTALLINE  ICE  <49p.) 

-  A  RATE-SENSITIVE  CONTINUUM  DAMAGE  MODEL  FOR  ICE 
(45p.)  -  Draft  Manuscript 


1 


3 


3.  TECHNICAL  REPORT 

A.  STATEMENT  OF  PROBLEM 

Introduction ;  Rate  mechanisms  govern  the  mechanical  behavior 
of  construction  materials  under  three  conditions:  \l)  at  high 
homologous  temperatures;  '(2)  when  subjected  to  short  term  loading 
that  is  quasi-static,  impulsive,  or  vibratory;  and  ^3)  when 
subjected  to  sustained  (or  creep)  loading  over  long  times. 
Examples  of  advanced  construction  materials  satisfying  one  or 
more  of  these  conditions  include: 

(a)  Polymer  matrix  structural  composites  such  as  FRPs/ 

j 

(b)  Single-ply  roofing  membranes  made  of  elastomers^ 

^c)  Polycrystalline  ice  in  cold  regions  engineering^'  ^  ^ 
.(d)  Cementitious  composites  such  as  high-strength 

concretes  and  FRCs. 

The  deformation  and  progressive  failure  of  such  rate  sensitive 
materials  is  governed  by  three  primary  mechanisms:  flow, 
distributed  cracking,  and  localized  crack  propagation.  The 
interactioi.  of  these  mechanisms  gives  rise  to  complex  mechanical 
behavior  on  the  macro-scale.  The  characterization  of  these 
mechanisms  and  their  interaction  is  therefore  of  fundamental 
importance.  ■ - 7  '^'-'27 

Any  given  material  may  display  purely  ductile,  purely 
brittle  or  combined  behavior  depending  upon  the  temperature  and 
conditions  of  loading.  The  primary  mechanism  of  flow  and  the 
constitutive  framework  of  rate  process  theory  are  appropriate  for 
characterizing  purely  ductile  behavior.  The  primary  mechanism 
associated  with  distributed  cracking  and  the  constitutive 
framework  of  damage  theory  are  appropriate  for  characterizing 
deformations  during  ductile-to-brittle  transition  or  when  the 
material  is  purely  brittle.  Such  deformations  are  accompanied  by 
the  formation  and  stable  growth  of  cracks  and/or  voids.  The 
primary  mechanism  associated  with  localized  cracking  and  fracture 
mechanics  theory  are  appropriate  for  predicting  the  onset  of 
material  instability. 

Damage  processes  control  failure  in  engineering  materials 


prior  to  localization  of  fracture  (formation  and  propagation  of  a 
a  single  crack).  Consequently,  the  ultimate  strength  of  such 
materials  is  generally  governed  by  damage.  In  addition,  damage 
plays  an  important  role  in  defining  the  onset  of  material 
instability  during  localized  fracture.  The  phenomenon  of  damage 
is  particularly  significant  when  the  state  of  stress  involves 
compression . 

Objectives  of  Research:  The  objective  of  this  research  is  to 
understand  and  characterize  rate  and  damage  processes  and  their 
interaction  in  advanced  construction  materials.  Both  theoretical 
and  experimentdl  research  are  being  pursued.  In  particular,  the 
following  three  tasks  are  being  undertaken: 

(a)  The  modeling  of  both  transient  and  steady  state  flow  in 
orthotropic  materials  based  on  physical  processes  occurring  on 
the  micro-scale. 

(b)  The  development  of  an  anisotropic  damage  model  which 
takes  into  account  the  damage  processes  leading  to  progressive 
deterioration  of  material  integrity.  Specifically,  the  rate  and 
pressure  sensitivity  of  the  damage  processes  are  investigated. 

(c)  The  development  of  an  experimental  capability  for 
testing  materials  at  low  temperatures,  under  high  rates  of 
loading,  and  under  high  confining  pressures.  The  program  of 
experimentation  is  being  designed  to  verify  the  theoretical 
models  of  flow  and  damage. 

Three  classes  of  materials  are  under  consideration: 
homogeneous  and  isotropic;  homogeneous  but  anisotropic; 
non-homogeneous ,  but  isotropic  (e.g.,  random  short-fiber 
reinforced  composites). 

B.  SUMMARY  OF  RESEARCH  ACCOMPLISHMENTS 

Research  during  the  past  year  has  led  to  the  development  of 
(i)  a  fflultiaxial  differential  flow  model,  (ii)  a  rate  and 
pressure  sensitive  anisotropic  damage  model,  and  (iii) 
specifications  for  testing  materials  as  described  Section  4. A. 

The  technical  approach  and  the  most  important  results  are 
summarized  here.  The  details  are  contained  in  Appendices  A  and  B 
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which  are  technical  papers  written  for  the  flow  and  damage 
models,  respectively,  with  specific  applications  to 
polycrystalline  ice.  This  material  has  been  chosen  for  the 
initial  phase  of  this  research  since  it  displays  a  wide  spectrum 
of  behaviors  typical  of  many  other  materials.  The  dominant 
behavioral  mode  is  strongly  dependent  on  the  rate  and  history  of 
loading,  as  well  as  the  ambient  temperature.  Furthermore, 
isotropic  polycrystalline  ice  is  an  ideal  testing  material  for 
verifying  the  damage  model  due  to  its  (i)  transparency  which 
facilitates  visual  correlation  of  cracking  activity  with 
quantitative  acoustic  emission  studies,  and  (ii)  homogeneity  and 
isotropy  which  lend  itself  easily  to  the  basic  verification  of 
theoretical  aspects  such  as  damage  induced  anisotropy  predicted 
by  the  evolution  laws  of  the  model. 

MULTIAXIAL  DIFFERENTIAL  FLOW  MODEL 

Physical  Basis. —  Flow  (or  creep)  is  modeled  in  terms  of  the 
physical  mechanisms  operating  on  the  micro-scale.  In 
polycrystalline  ice,  it  is  known  from  both  theoretical  and 
experimental  work  that  at  least  two  deformation  systems,  a  soft 
system  and  a  hard  system,  are  present  during  flow.  These  include 
grain  boundary  sliding  and  basal  slip,  or  basal  slip  and  slip  on 
a  non-basal  plane. 

Under  load  the  material  initially  resists  the  applied  stress 
in  an  elastic  manner  and  then  flow  begins  on  the  soft  and  hard 
systems.  Flow,  particularly  on  the  easy  soft  system,  causes  the 
build  up  of  internal  elastic  stresses.  This  may  occur  as  a  result 
of  grain  boundary  sliding  next  to  grains  poorly  aligned  for 
deformation  or  dislocation  pile-ups  at  the  boundaries  of  such 
grains.  These  internal  elastic  stresses  are  known  as  the  back  cr 
rest  stresses.  In  addition,  internal  drag  stresses  which  resist 
dislocation  fluxes  and  sliding  are  generated  by  (a)  creep 
resistant  substructures,  i.e.  subgrains  and  cells,  formed  by 
grain  boundary  sliding  or  dislocation  movement,  and  (b) 
dislocation  entanglement,  dipole  formation,  and  kink  band 
formation  during  slip. 


6 


An  increasing  back  stress  contributes  to  kinematic  hardening 
which  induces  directionally  dependent  material  behavior,  referred 
to  as  deformation  induced  anisotropy.  The  Baushinger  effect  in 
metals  is  an  example  of  kinematic  hardening.  On  the  other  hand, 
an  increasing  drag  stress  con^ributes  to  isotropic  hardening.  In 
isotropic  hardening,  subsequent  material  properties  are 
independent  of  the  direction  of  pre-straining. 

In  the  model  the  total  deformation  resulting  from  the 
interaction  between  the  soft  and  hard  systems  is  decomposed  into 
a  transient  component  and  a  steady  state  component.  Isotropic  and 
kinematic  hardening  phenomena  are  active  during  transient  flow 
and  give  rise  to  elastic  strains  which  are  recoverable  upon 
unloading.  These  time-dependent  elastic  strains  represent  delayed 
elasticity  or  anelasticity .  Steady  state  flow  is  associated  with 
viscous  strains  which  are  irrecoverable.  Such  viscous  deformation 
is  attributed  to  intragranular  deformation  processes,  especially 
to  the  movement  of  dislocations.  When  transient  flow  has  become 
saturated,  hardening  effectively  ceases  and  pure  viscous  flow 
results . 

Mathematical  Formulation. —  All  equations  stated  in 
connection  with  the  Multiaxial  Differential  Flow  Model  refer  to 
those  in  Appendix  A.  The  uniaxial  governing  equation  is  obtained 
by  expressing  the  total  strain  rate  as  a  sum  of  the  transient  and 
steady  state  flow  components,  as  well  as  the  elastic  component 
associated  with  instantaneous  deformation  (Eq.  (1)).  The 
instantaneous  component  is  described  by  the  classical  theory  of 
linear  elasticity  which  in  rate  form  relates  the  elastic  strain 
rate  to  the  stress  rate  through  the  Young's  modulus  (Eq.  (2)). 

The  viscous  strain  follows  the  well-known  power  law  of  Glen 
(1955),  i.e.  the  viscous  strain  rate  is  proportional  to  the  third 
power  of  stress  (Eq.  (3)).  The  constant  in  Glen's  law  is 
temperature  dependent  and  follows  the  Arrhenius  activation  energy 
law  (Eq.  (4)).  The  transient  strain  rate  is  also  taken  to  follow 
Glen's  law  with  the  same  temperature  dependence,  but  transient 
flow  is  now  driven  by  a  reduced  stress  quantity,  equal  to  the 
ratio  of  the  applied  stress  less  the  back  stress  to  a 
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dimensionless  drag  stress  (Eq.  (5)). 

Since  the  back  and  drag  stresses  are  history  dependent 
variables  associated  with  transient  flow,  evolution  equations 
must  be  specified  for  them.  Due  to  the  elastic  nature  of  the 
recoverable  transient  strains,  the  time  rates  of  change  of  the 
back  and  drag  stresses  can  be  assumed  to  be  proportional  to  the 
transient  strain  rate  itself.  These  evolution  equations  are 
stated  in  Eqs.  (6)  and  (7). 

Under  creep  loading  Eq.  (6)  predicts  that  the  back  stress 
will  asymptotically  increase  to  a  value  equal  to  the  applied 
stress,  at  which  point  transient  flow  will  cease.  Under  constant 
strain  rate  loading  the  back  stress  will  asymptotically  attain 
the  value  of  the  steady-state  stress.  For  either  loading 
condition,  the  drag  stress  will  also  reach  its  limiting  value 
when  transient  flow  ceases.  Under  reversed  or  cyclic  loading  the 
back  stress  will  reverse  or  switch  back  between  positive  and 
negative  values,  i.e.,  the  physical  processes  associated  with 
kinematic  hardening  can  locally  relax  or  move  back  and  forth 
between  positive  and  negative  values.  During  reversed  loading, 
the  drag  stress  will  also  decrease,  indicating  a  decreasing 
resistance  to  transient  flow.  On  the  other  hand,  the  drag  stress 
will  reach  a  saturated  value  under  cyclic  loading,  while  the  back 
stress  will  eventually  switch  back  and  forth  between  some 
limiting  values. 

Model  Formulation  in  Dimensionless  Variables. —  Based  on  the 
test  data  of  Jacka  (1984)  for  isotropic  polycrystalline  ice, 

Ashby  and  Duval  (1985)  have  suggested  that  unique  relationships 
exist  between  certain  dimensionless  variables.  Such  relationships 
are  predicted  by  the  proposed  model  as  explained  below. 

Under  constant  stress  loading  the  model  predicts  that  the 
relationships  between  dimensionless  strain,  strain  rate  and  time 
are  independent  of  the  applied  stress  and  temperature.  The 
dimensionalization  is  achieved  by  expressing  the  dimensionless 
variables  as  ratios  between  the  variables  and  their  appropriate 
normalizing  variables,  i.e.,  all  strain  quantities  are  normalized 
with  the  elastic  strain,  all  strain  rate  quantities  with  the 
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steady  viscous  strain  rate,  and  all  stress  quantities  with  the 
applied  stress  (Eqs.  (8),  (9),  (11)),  Dimensionless  time  is 
defined  as  time  multiplied  by  the  viscous  strain  rate  divided  by 
the  elastic  strain  ( Eq .  (10)).  It  is  possible  to  express  the 
equations  of  the  model  as  well  as  the  evolution  equations  in 
terms  of  these  dimensionless  variables  (Eqs.  (12)-(17)),  thus 
predicting  unique  relationships  (i.e.,  independent  of  the  applied 
stress  and  temperature)  between  them. 

In  the  case  of  constant  strain  rate  loading,  the  model 
predicts  that  a  unique  relationship  exists  between  the 
dimensionless  stress  and  time  (or  strain).  The  dimensi onal i za t i on 
is  achieved  by  normalizing  all  strain  variables  with  the  elastic 
strain,  and  all  stress  variables  with  the  steady-state  stress 
(Eqs.  (18),  (20)).  The  steady-state  stress  is  given  by  Glen's 
power  law  for  viscous  flow,  Eq.  (21).  Dimensionless  time  is 
defined  as  time  multiplied  by  the  constant  strain  rate  divided  by 
the  maximum  elastic  strain  ('^.q.  (19)).  The  equations  of  the  model 
and  the  governing  equations  can  be  expressed  in  terms  of  these 
dimensionless  variables  (Eqs.  (22)-(27)).  Since  the  temperature 
dependent  constant  associated  with  Glen's  law  and  the  applied 
constant  strain  rate  have  been  eliminated  from  these  equations,  a 
unique  master  curve  can  be  used  to  relate  the  dimensionless 
stress  and  time. 

Experimental  Validation  of  Uniaxial  Model. —  The  uniaxial 
model  is  verified  against  the  constant  stress  (creep)  data  of 
Jacka  (1984)  and  Sinha  (1978).  The  former  data  set  is  obtained 
from  tests  conducted  on  isotropic  polycrystalline  ice  with  an 
average  grain  size  of  1.7  mm,  whereas  the  latter  is  obtained  from 
tests  conducted  on  columnar-grained  ice  with  an  average  grain 
size  of  3  BUB.  To  further  demonstrate  the  capability  of  the  model, 
the  predicted  strain  response  under  a  monotonically  increasing 
stress  is  verified  against  Sinha's  (1981)  data  for 
columnar-grained  ice  with  an  average  grain  size  of  4.5  mm. 

Finally,  the  creep  and  recovery  response  of  randomly  oriented 
snow  ice  is  studied  using  Brill  and  Camp's  data  taken  from  Sinha 
(1979).  Note  that  all  figures  mentioned  in  this  subsection  refer 
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to  those  in  Appendix  A. 

Figures  1-3  show  Jacka's  data  together  with  the  model 
predictions  (curves  (b))  using  dimensionless  variables.  As  can  be 
seen  from  these  figures,  unique  relationships  between 
dimensionless  strain  rate  and  time,  between  dimensionless  strain 
rate  and  strain,  and  between  dimensionless  strain  and  time,  are 
fully  justified  by  the  test  data  corresponding  to  various  stress 
levels  and  temperatures. 

Excellent  agreement  between  model  prediction  and  Sinha's 
creep  data  for  various  temperatures  is  demonstrated  in  Fig.  4. 
Under  a  monotonically  increasing  stress  history,  the  predicted 
strain-time  and  stress-strain  relationships  also  agree  well  with 
Sinha's  test  data,  as  shown  in  Fig.  6. 

Figure  7  shows  three  sets  of  creep  and  recovery  data 
obtained  under  different  conditions  of  stress  and  temperature. 

The  model  predictions  for  both  creep  and  recovery  agree  fairly 
well  with  the  data. 

Multiaxial  Model  Formulation. —  The  three-dimensional 
generalization  of  the  uniaxial  model  follows  from  the  uniaxial 
formulation,  i.e.,  it  is  based  on  strain  decomposition,  linear 
elasticiy,  and  the  rate  theory  of  flow.  Constitutive  relations 
are  derived  for  each  mechanism  of  deformation  in  the  model, 
resulting  in  the  orthotropic  equivalent  of  Eqs.  (l)-(7). 

Orthotropic  elasticity  (Eq.  (44))  follows  the  classical 
theory  (see,  for  example,  the  book  by  Lekhnitskii,  1963).  The 
orthotropic  generalization  of  viscous  flow  is  derived  from  the 
normality  of  the  viscous  strain  rate  to  a  scalar-valued  viscous 
flow  potential  defined  in  terms  of  an  equivalent  stress  measure 
(see  Eqs.  (45)-(52)).  The  equivalent  stress  measure  with  six 
orthotropic  constants  is  a  generalized  form  for  orthotropic 
materials.  An  equivalent  strain  or  strain  rate  measure  (Eq.  (56)) 
for  orthotropic  materials  can  also  be  derived  based  on  the 
hypothesis  of  energy  equivalence.  The  orthotropic  genral i zation 
of  transient  flow  is  based  on  the  assumption  of  flow 
incompressibility.  Furthermore,  the  material  orthotropy  is 
described  by  the  same  set  of  orthotropic  constants  mentioned 
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previously.  In  a  manner  similar  to  the  derivation  of  the  viscous 
strain  rate,  the  equation  for  transient  flow  is  derived  from  the 
normality  of  the  transient  strain  rate  to  a  scalar-valued 
transient  flew  potential  defined  in  terms  of  an  equivalent 
reduced  stress  measure.  Note  that  since  incompressibility  of  flow 
is  assumed,  all  str  s  and  strain  quantities  used  in  the 
equations  for  viscous  and  tranient  deformations  are  deviatoric. 
The  resulting  expression  is  expressed  in  Eq.  (58).  Evolution 
equations  are  also  generalized  to  three  dimensions  where  the 
apppropriate  variables  are  the  deviatoric  bacit  stress  vetor  and  a 
scalar  equivalent  drag  stress  (Eqs.  (63)  and  (64)).  The 
constraint  condition  (i.e.,  strain  decomposition),  the  equations 
for  orthotropic  elasticity,  viscous  and  transient  strain  rates, 
together  with  the  evolution  equations  form  the  equations  of  the 
multiaxial  model  which  can  be  numerically  integrated  to  predict 
the  model  response  under  variable  loading  histories. 

Experimental  Validation  of  Multiaxial  Model. —  In  this 
subsection  the  x-axis  is  taken  to  be  normal  to  the  columnar  ice 
sheet  which  is  defined  by  the  y-z  plane.  The  columns  of  the  ice 
sheet  are  parallel  to  the  growth  direction,  i.e.  along  the 
x-axis.  The  c-axes  of  the  ice  crystals  are  assumed  to  lie  in  the 
y-z  plane  and  are  aligned  in  the  y-direction. 

The  model  response  is  verified  against  Frederking's  (1977) 
data  obtained  from  plane  strain  compression  tests  on 
columnar-grained  transversely  isotropic  (a  special  case  of 
orthotropy  in  which  parallel  planes  of  isotropy  exist  in  the 
material)  freshwater  ice.  The  planes  of  isotropy  are  the  y-z 
planes.  The  orthotropic  constants  can  be  determined  from  uniaxial 
tests  (see  Eq.  (69)  for  details).  In  his  type  A  tests,  strains  in 
the  z-direction  are  constrained  to  zero  and  stresses  are  applied 
in  the  y-direction.  At  steady  state  where  the  power  law 
orthotropic  formulation  suffices,  the  ratio  of  the  plane  strain 
stress  to  the  unconfined  stress  at  the  same  strain  rate  (Eq.  74) 
is  predicted  to  vary  between  2. 1-5.1.  These  values  are  consistent 
with  Frederking's  experimental  observations  which  were  close  to  2 
at  high  strain  raes  and  to  5  at  low  strain  rates.  In  his  type  B 
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tests,  strains  in  the  x-direction  are  constrained  to  zero  while 
stresses  are  again  applied  in  the  y-direction.  The  predicted 
ratio  between  the  plane  strain  stress  and  the  unconfined  stress 
at  the  same  strain  rate  (Eq.  75)  varies  between  1.01  and  1.06, 
which  again  agrees  with  Frederking's  test  results  which  showed 
negligible  influence  of  x-firection  confinement  on  stresses. 

RATE-SENSITIVE  CONTINUUM  DAMAGE  MODEL 

Rationale  For  the  Use  of  Damage  Variable  in  Constitutive 
Modeling. —  The  theory  of  plasticity,  originally  developed  as  a 
phenomenological  theory  for  modeling  nonlinear  material  behavior 
due  to  such  irreversible  processes  as  crystalline  slip  and 
twinning,  is  generally  found  to  be  unsuitable  for  the  description 
of  material  behavior  dominated  by  the  nucleation  and  growth  of 
microcracks.  Conventional  fracture  mechanics,  on  the  other  hand, 
is  concerned  with  the  prediction  of  material  response  dominated 
by  a  single  macrocrack.  However,  there  are  many  instances  when 
the  material  deteriorates  by  the  nucleation  and  stable  growth  of 
a  multitude  of  stable  microcracks  (damage)  distributed  over  the 
entire  material  continuum,  especially  when  the  material  is  loaded 
in  compression.  Depending  on  the  rate  of  loading  and  the  stress 
state,  the  evolution  of  damage  may  significantly  affect  material 
behavior  and  the  flow  model  by  itself  is  therefore  insufficient 
for  a  proper  description  of  material  behavior. 

The  prediction  of  the  development  of  damage  under  different 
strain  rates  and  stress  states  is  thus  of  considerable  importance 
if  an  accurate  prediction  of  material  behavior  is  to  be  obtained. 
Since  it  is  desired  to  study  the  development  of  damage  per  se, 
isotropic  polycrystalline  ice  is  chosen  as  a  material  for  study 
where  the  complexity  of  material  anisotropy  is  absent. 
Furthermore,  ice  also  displays  rate  and  pressure  sensitive 
behavior . 

During  the  creep  damage  process,  the  nucleation  and  growth 
of  microcracks  can  be  influenced  by  diffusion,  by  dislocation 
creep,  and  by  stress  concentrations  at  irregularities  on  grain 
boundaries  caused  by  grain  boundary  sliding  (Sinha,  1984).  These 
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processes  are  linked  to  the  development  of  creep  strain,  which 
therefore  affects  damage  growth.  On  the  other  hand,  damage  can  be 
suppressed  by  a  confining  pressure.  Under  low  strain  rate 
loading,  microcraching  is  insignificant  and  the  application  of  a 
hydrostatic  state  of  confining  stress  does  not  increase  the 
maximum  shear  stress  that  can  be  sustained  by  the  material.  Under 
high  strain  rate  loading,  damage  is  significant  and  application 
of  pressure  suppresses  the  damage,  as  a  result  of  which  a  higher 
shear  strength  can  be  sustained  by  the  material.  The  ratio  of  the 
confined  shear  strength  to  the  unconfined  shear  strength  in  the 
case  when  there  is  considerable  damage  supression  is  greater  than 
that  in  the  case  when  almost  pure  flow  with  little  damage 
dominates  the  material  response.  This  phenomenon  is  known  as 
pressure  sensitivity.  Under  very  large  confining  pressures  the 
material  may  undergo  pressure  melting,  and  its  ability  to  sustain 
shear  decreases  and  eventually  vanishes  at  the  pressure 
corresponding  to  phase  change. 

Damage  in  ice  leads  to  "strain  softening"  under  constant 
strain  rate  loading,  and  to  accelerating  or  tertiary  creep  under 
constant  stress  loading.  Since  a  reduction  in  elastic  modulus  is 
observed  during  unloading  (Uordaan,  1986),  it  can  be  inferred 
that  damage  affects  the  elastic  (in  addition  to  creep)  properties 
'jf  ice  as  well.  Furthermore,  the  microcracks  that  have  been 
observed  in  ice  can  be  approximated  as  planar  in  shape,  and  they 
are  oriented  with  their  normals  in  the  positive  principal  stress 
directions  (Sinha,  1982).  Such  induced  anisotropy  is  also 
observed  in  many  materials  and  is  especially  important  under 
non-proportional  loading  when  conditions  for  rupture  or  failure 
depends  on  the  particular  loading  history  of  the  material,  and 
hence  the  orientation  of  the  microcracks  (Chaboche,  1982, 

Murakami  and  Ohno,  1981).  Consequently,  in  modeling  damage  a 
scalar  representation  is  inadequate  and  a  vectorial  or  tensorial 
approach  is  necessary. 

The  model  presented  here  is  based  on  the  assumption  of  an 
initially  isotropic  material  and  takes  into  account  strain  rate 
and  pressure  sensitivity,  pressure  melting  and  damage-induced 
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anistropy.  It  uses  a  second  order  damage  tensor  representation 
introduced  by  Murakami  and  Ohno  (1981).  Damage  effects  are 
incorporated  through  the  use  of  a  net  stress  tensor.  This  is 
based  on  the  concept  of  locally  magnified  stresses  arising  from 
a  reduction  in  load  bearing  area  due  to  damage.  The  evolution 
equations  consider  damage  induced  by  both  local  tensile  and  shear 
stresses.  The  damage  model  is  incorporated  in  the  flow  model 
through  the  use  of  an  averaged  net  stress  tensor,  i.e.,  the 
stress  tensor  in  the  flow  model  is  replaced  by  the  averaged  net 
stress  tensor. 

Continuum  Damage  Modeling. —  To  predict  material  response 
when  damage  effects  are  not  negligible,  it  is  necessary  to  (i) 
select  an  appropriate  damage  variable  based  on  some  averaging 
over  a  representative  volume,  (ii)  establish  the  damage  evolution 
laws,  and  (iii)  formulate  an  appropriate  set  of  constitutive 
relations  which  incorporate  the  effects  of  damage.  Note  that  all 
equations  below  refer  to  those  in  Appendix  B. 

The  microcracks  are  directly  used  to  define  the  damage 
tensor  through  the  reduction  in  the  strength  of  sections.  This  is 
therefore  a  geometric  approach.  For  an  elementary  volume  of  the 
damaged  material  continuum  associated  with  a  material  point,  it 
has  been  shown  that  the  damage  state  can  be  defined  by  a  second 
order  symmetric  tensor  (Eq.  (1)).  This  damage  tensor  is 
constructed  by  taking  the  tensor  product  of  the  k^^  unit  vector 
normal  to  the  k^^  planar  microcrack  with  itself,  integrating  over 
the  area  of  the  microcrack,  and  summed  for  all  microcracks  in  the 
elementary  volume.  Since  the  resulting  tensor  is  second  order 
symmetric,  it  can  be  written  in  terms  of  its  three  principal 
values  and  principal  directions  (Eq.  (2)),  i.e.,  as  a  sum  of  the 
three  tensor  products  between  the  unit  vectors  representing  the 
principal  directions  with  themselves,  where  each  tensor  product 
is  multiplied  by  the  corresponding  principal  value.  This  implies 
that  an  arbitrary  damage  field  can  be  described  in  terms  of  three 
mutually  orthogonal  systems  of  parallel  microcracks.  The  effect 
of  a  multitude  of  damage  fields  is  in  general  obtained  by  summing 
the  different  damage  tensors  (one  for  each  field). 
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A  net  stress  tensor  can  be  used  to  represent  the  locally 
magnified  stress  which  governs  the  evolution  of  damage. 
Mathematically  it  can  be  defined  as  the  symmetric  part  of  the 
Cauchy  stress  multiplied  by  a  damage  effect  tensor  (see  Eq . ( 8 ) ) . 
The  damage  effect  tensor  is  equal  to  the  inverse  of  the 
difference  beween  the  identity  tensor  and  the  damage  tensor  (see 
Eq.  (9).)  Physically  the  damage  effect  tensor  can  be  thought  of 
as  a  quantity  representing  the  effect  of  local  stress 
concentration  which  influences  damage  growth.  Since  at  the 
present  the  emphasis  is  on  the  study  of  how  damage  evolves  and 
not  how  damage  affects  the  three  dimensional  deformation 
behaviour,  an  averaged  net  stress  is  used  to  specify  the 
constitutive  relations.  This  averaged  tensor  is  defined  as  the 
mean  of  the  diagonal  components  of  the  damage  effect  tensor 
multplied  by  the  stress  tensor  (Eq.  (10)). 

Evolution  Equations  of  Damage. —  The  net  stress  tenser  is 
used  to  specify  the  evolution  laws.  To  avoid  cumbersome 
terminology,  it  is  understood  that  in  this  subsection  all 
stresses  refer  to  net  stresses.  It  has  been  experimentally 
observed  for  many  materials  that  microcracks  tend  to  form  in 
planes  perpendicular  to  the  local  principal  tensile  stresses 
(Murakami  and  Ohno,  1981;  Costin,  1983).  Under  a  non-hydrostatic 
compressive  state  of  stress  (which  is  the  stress  state  considered 
here),  local  tensile  stresses  can  arise  fom  material  property 
mismatch  between  grains  or  from  contact  stresses  between  grains 
with  irregular  grain  boundaries  (Costin,  1983).  Since  the  normals 
to  the  microcracks  formed  under  compression  are  usually  oriented 
in  the  directions  of  the  principal  tensile  deviatoric  stresses, 
it  can  be  assumed  that  the  local  tensile  stresses  act  in  the  same 
directions.  In  the  proposed  model  it  is  further  assumed  that  such 
microcracks  form  only  in  the  planes  perpendicular  to  the  maximum 
local  principal  tensile  deviatoric  stress.  The  "directionality" 
of  anisotropic  damage  associated  with  microcracks  formed  under 
local  tension  can  then  be  determined  by  the  tensor  product  of  the 
unit  vector  in  the  direction  of  the  maximum  local  principal 
tensile  deviatoric  stress  with  itself. 
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The  spherical  or  hydrostatic  state  of  stress  is  known  to 
suppress  the  formation  and  growth  of  microcracks.  This  phenomenon 
is  also  observed  in  ice  (Mellor,  1983),  especially  at  high  rates 
of  loading.  In  the  evolution  laws  this  effect  is  modeled  by 
reducing  the  damage  due  to  the  maximum  local  stress  by  an  amount 
proportional  to  the  hydrostatic  stress. 

In  addition  to  the  maximum  local  tensile  stress,  local  shear 
stresses  may  also  generate  microcracks.  Such  local  stresses  may 
be  regarded  as  related  to  the  second  invariant  of  the  deviatoric 
stress  tensor,  and  the  microcracks  generated  in  this  manner  give 
rise  to  off-diagonal  components  in  the  dama'^e  tensor.  However, 
the  damage  tensor  can  be  transformed  to  give  principal  normal 
values  on  the  main  diagonal  only.  The  directionality  associated 
with  the  second  invariant  of  the  stress  deviator  may  be 
approximately  regarded  as  isotropic  since  the  majority  of  the 
shear  cracks  tend  to  be  randomly  orientated  with  respect  to  the 
global  frame  and  will  consequently  give  rise  to  equal  components 
in  all  three  principal  directions.  Thus,  local  shear  stresses  are 
assumed  to  contribute  to  ’’isotropic"  damage. 

To  take  into  account  the  effect  of  strain  rate  on  damage,  it 
is  proposed  that  a  power  law  holds  between  damage  growth  rate  and 
the  effective  total  strain  rate  (Murakami  et  al.,  1986). 
Furthermore,  the  effects  of  local  stress  concentration  on  damage 
evolution  are  modeled  by  using  a  scalar  invariant  of  the  damage 
effect  tensor  for  the  isotropic  part  of  the  damage  evolution 
laws,  and  by  using  the  magnitude  of  the  damage  effect  in  the 
direction  of  the  maximum  principal  stress  deviator  for  the 
anisotropic  part.  The  complete  evolution  equations  are  given  by 
Eqs.  (13)-(15). 

Constitutive  Relations. —  The  constitutive  laws  used  here 
are  based  on  a  simple  nonlinear  generalization  of  the  Maxwell 
model.  Details  are  given  in  Appendix  B,  where  the  constitutive 
relations  are  given  by  Eqs.  (19)-(22).  The  flow  model  described 
earlier  is  in  the  process  of  being  implemented  and  will  supercede 
the  Maxwell  model.  The  essential  concept  in  incorporating  damage 
effects  in  the  constitutive  relations,  however,  is  to  replace  the 
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stress  components  by  appropriate  net  stress  components.  As 
described  earlier,  averaged  net  stress  components  are  used, 
implying  that  the  average  damage  effect  on  the  global  stress  is 
isotropic.  This  is  only  true  when  damage  effect  is  small 
(Murakami  and  Imaizumi,  1982).  Research  is  being  undertaken  to 
use  a  better  representation  of  the  form  of  the  net  stress  for  use 
in  the  constitutive  relations. 

Pressure  Melting  Model. —  Ice  is  known  to  change  phase  at  a 
certain  level  of  hydrostatic  stress.  Triaxial  testing  of  ice  was 
carried  out  relatively  recently  (see,  for  example,  Jones,  1982). 
It  is  experimentally  observed  that  the  resistance  of  ice  to  shear 
stress  decreases  when  hydrostatic  stress  reaches  some  moderate 
level.  Thus  the  flow  resistance  parameter  in  Glen's  power  law 
should  decrease  with  increase  in  hydrostatic  stress.  On  the  other 
hand,  it  is  also  known  that  this  parameter  will  decrease  with 
increase  in  temperature  according  to  the  Arrhenius  law.  Thus,  its 
variation  with  pressure  can  be  modeled  through  a  temperature 
correction  using  the  phase  diagram  for  ice.  The  phase  diagram 
plots  the  pressure-temperature  relationship  corresponding  to 
phase  change.  The  relation  between  the  flow  parameter  and  the 
corrected  temperature  is  then  obtained  using  the  Arrhenius  law. 

Experimental  Validation. —  The  damage  model  is  verified 

against  several  independent  sets  of  data,  including  the  uniaxial 

data  of  Gold  reproduced  in  Mellor  (1983),  and  the  triaxial  data 

of  Jones  (1982).  The  maximum  stress  is  plotted  against  the  strain 

rate  in  Fig.  7,  which  shows  Gold's  data  for  various  ice  types 

tested  under  uniaxial  compression  at  -10”C,  as  well  as  the  model 

prediction.  The  agreement  between  theory  and  data  is  excellent. 

These  results  show  that  damage  is  insignificant  at  strain  rates 
—4  —1 

below  about  10  s  where  deformation  is  dominated  by  pure  flow. 
Thus  the  maximum  stress  or  strength  of  the  material  increases 
with  strain  rate  according  to  Glen's  power  law.  For  strain  rates 
greater  than  10”^  s”^ ,  i.e.,  in  the  ductile-to-bri ttle  transition 
region,  damage  becomes  significant  and  the  strength  becomes  less 
than  that  predicted  by  the  power  law,  but  still  increases 
monotonically  until  a  plateau  is  reached  at  very  high  strain 


17 


rates . 

The  data  of  Jones  for  isotropic  polycrystalline  ice  is 
obtained  from  triaxial  tests  conducted  at  -12®C.  The  confining 
pressure  and  the  strain  rate  of  compressive  loading  are  varied 
in  order  to  study  their  effects  on  the  shear  strength.  The  data 
and  the  model  predictions  agree  well,  thus  validating  the  damage 
model  with  regard  to  the  effects  of  strain  rate,  stress  state  and 
pressure  melting.  Specifically,  these  results  show  that  (i)  shear 
strength  increases  with  the  applied  constant  strain  rate,  (ii) 
shear  strength  increases  with  moderate  confining  pressure,  and 
(iii)  pressure  melting  lowers  the  shear  resistance.  Furthermore, 
the  damage  theory  predicts  that  damage  is  predominantly 
anisotropic  (i.e.,  microcracks  are  preferentially  oriented 
perpendicular  to  the  maximum  principal  stress  deviator)  when 
there  is  no  confinement,  but  a  moderate  level  of  confining  stress 
causes  the  formation  of  shear  microcracks,  giving  isotropic 
damage.  Anisotropic  damage  becomes  relatively  less  significant 
with  increase  in  confining  pressure.  Ultimately  at  very  high 
pressure,  damage  growth  is  almost  entirely  associated  with  the 
second  invariant  of  the  stress  deviator,  and  the  material 
undergoes  pseudo-flow  with  many  distributed  microcracks. 

EXPERIMENTAL  FACILITIES  AND  PROGRAM 

The  experimental  facilities  includ'e  a  low  temperature 
materials  testing  facility,  a  load  frame  and  a  triaxial  cell 
required  for  mechanical  testing,  acoustic  emission  monitoring 
system  with  host  computer,  transducers  and  multi-channel 
monitoring  system,  and  a  data  acquisition  system  composed  of  a 
microprocessor/controller  and  a  high  speed  analog  to  digital 
converter.  Part  of  the  equipment  has  been  purchased  and  installed 
with  partial  ARO  support  during  1987.  Additional  funding  has  been 
requested  from  ARO  for  the  remaining  equipment  needs  (load  frame, 
triaxial  cell,  acoustic  emission  monitoring  system.) 

Low  Temperature  Materials  Testing  Facility. —  This  consists 
of  three  cold  rooms,  a  test  chamber,  a  growth  room  and  a  storage 
room.  The  minimum  temperature  in  the  test  chamber  can  reach  -40®C 
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and  the  precision  of  the  temperature  can  be  controlled  to  within 
+1®C.  The  growth  room  is  used  for  growing  ice  samples  by  freezing 
distilled,  degassed  and  deionized  water  at  0*C  through  seed  ice 
crystals.  The  temperature  range  in  this  room  lies  between  -10  and 
0®C  and  the  temperature  precision  is  within  +1®C.  The  storage 
room  is  used  for  storing  ice  samples  and  seed  crystals,  and  as  a 
buffer  between  the  other  two  rooms  and  the  exterior.  Its 
temperature  range  is  the  same  as  that  for  the  growth  room,  but 
the  temperature  precision  is  coarser,  within  +2®C. 

Load  Frame  and  Triaxial  Cell. —  A  load  frame  capable  of 
delivering  220  kips  at  a  displacement  rate  of  1  in/sec  has  been 
specified.  The  frame  can  also  have  vertical  stress  control.  The 
triaxial  cell  will  be  designed  for  a  maximum  confining  pressure 
of  11000  psi.  The  confining  stress  will  be  applied  in  fixed 
proportion  to  the  vertical  axial  stress.  This  requires  strain 
control  of  axial  loading  with  a  second  loop  to  adjust  the 
confining  pressure. 

Acoustic  Emission  Monitoring  System. —  The  system  will  be 
capable  of  performing  three-dimensional  location  analysis  of 
microcracks  using  six  resonant  transducers  located  around  a 
cylindrical  specimen.  Transducers  should  have  a  resonant 
frequency  of  about  0.2  MHz  and  a  frequency  bandwidth  of  between 
0.1  to  0.6  MHz.  Discrimination  time  is  at  most  150  micro-seconds. 
Fast  data  acquisition  together  with  the  ability  to  dump  data 
directly  onto  disk  will  allow  for  rapid  data  and  post  test 
analyses . 

The  experimental  program  has  been  designed  to  verify  the 
theoretical  aspects  of  the  flow  and  damage  models,  i.e.,  (i)  to 
generate  a  comprehensive  set  of  experimental  data  at  constant 
strain  rates  for  characterizing  the  deformation  of  materials  when 
damage  processes  are  active;  (ii)  to  investigate  the  formation  of 
"first"  cracks  during  deformation  by  monitoring  acoustic 
emissions;  (iii)  to  develop  and  apply  quantitative  acoustic 
emission  theory  for  locating  cracks  (position,  time,  direction, 
and  size);  and  (iv)  to  theoretically  characterize  the 
rate-sensitive  evolution  of  damage  during  deformation  by  flow. 
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All  four  areas  of  the  investigation  will  be  carr 
and  high  strain  rates  of  loading,  under  uniaxial 
loading  conditions,  and  at  various  temperatures, 
program  should  commence  towards  the  end  of  year 
1988)  . 
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1.  INTRODUCTION 

Boundary  value  problems  in  applied  ice  mechanics  involving 
multiaxial  states  of  stress  and  complex  loading  histories,  such 
as  those  encountered  during  ice-structure  interaction,  are 
increasingly  being  solved  using  numerical  models  including  the 
finite  element  method  (Jordaan,  1986).  Constitutive  models  are 
required  to  characterize  the  ice  deformation  by  viscoelastic  flow 
in  numerical  simulations. 

In  problems  where  only  "steady  state"  flow  is  of  interest, 
an  elastic  -  power  law  creep  model  of  ice  (sometimes  without  the 
elastic  component)  is  adequate.  The  most  widely  used  model  of 
steady  state  or  viscous  flow  cf  polycrystalline  ice  is  Glen's 
power  law.  The  multiaxial  generalization  of  the  differential 
model  follows  from  conventional  elasticity  theory  and  from  the 
rate  theory  of  flow.  The  latter  is  based  on  normality  of  the 
viscous  deformation-rate  to  a  scalar  valued  flow  potential 
expressed  in  terms  of  an  equivalent  stress  measure.  Palmer  (1967) 
has  derived  the  multiaxial  law  for  incompressible  flow  of 
isotropic  ice,  while  Shyam  Sunder,  Ganguly  and  Ting  (1987)  have 
presented  an  orthotropic  model  of  incompressible  flow. 

Both  the  elastic  and  "transient"  flow  behavior  of  ice, 
however,  are  of  great  importance  in  a  broad  range  of  ice 
mechanics  problems  (Gold,  1977,  Sinha  et  al.,  1987).  The  most 
widely  used  flow  law  for  ice  under  uniaxial  loading  is  the  creep 
compliance  function  proposed  by  Sinha  (1978,  1979).  This 
formulation  postulates  that  grain  boundary  sliding  governs 
transient  deformation,  and  that  the  complian'.e  function  is 
linearly  dependent  on  stress  and  nonlinearly  dependent  on  time. 
For  conditions  other  than  constant  stress  creep,  monotoni cally 
increasing  stress  in  particular,  Sinha  (1983)  has  applied  the 
nonlinear  compliance  f-nction  in  conjunction  with  a  convolution 
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integral  to  predict  the  mechanical  behavior.  This  integral 
formulation  assumes  a  particular  generalization  of  Boltzmann's 
superposition  principle  for  transient  deformations. 

Le  Gac  and  Duval  (1980)  have  proposed  multiaxial 
constitutive  relations  for  the  inelastic  deformation  of 
polycrystalline  ice  which  account  for  the  phenomena  of  isotropic 
and  kinematic  hardening.  Considering  deformation  mechanisms  in 
ice,  Ashby  and  Duval  (1985)  have  subsequently  developed  a 
kinematic  hardening  model  based  on  a  two-bar  truss  analogy.  They 
have  used  the  model  to  identify  certain  dimensionless  variables 
from  which  a  single  master  curve  can  be  developed  for  the  creep 
of  polycrystalline  ice.  The  appropriateness  of  the  variables  has 
been  demonstrated  using  the  comprehensive  experimental  data  of 
Jacka  (1984)  for  dense  isotropic  polycrystals  with  a  mean  grain 
size  of  1.7  mm.  The  predictive  capability  of  their  model, 
however,  has  not  been  explicitly  tested  against  this  data  set. 

This  paper  presents  a  differential  flow  model  for  the 
deformation  of  polycrystalline  ice  which  (i)  accounts  for  both 
isotropic  and  kinematic  hardening,  and  (ii)  satisfies  the 
dimensional  requirements  identified  by  Ashby  and  Duval  (1985). 
Flow  (or  creep)  is  modeled  in  terms  of  two  nonlinear 
def ormiation-rate  mechanisms;  the  first  mechanism  governs  the 
transient  deformation-rate  (creep)  which  decays  to  zero  as  both 
an  elastic  back  stress  and  a  drag  stress  measure  increase 
asymptotically;  the  second  mechanism,  which  is  modeled  in  terms 
of  the  well-known  power  law,  governs  the  viscous  deformation- 
rate.  The  evolution  of  the  back  or  rest  stress  contributes  to 
kinematic  hardening,  while  that  of  the  drag  stress  contributes  to 
isotropic  hardening. 

In  general,  numerical  integration  of  the  governing  equations 
is  necessary  for  predicting  the  model  response  under  arbitrary 
loading  histories  since  both  isotropic  and  kinematic  hardening 
are  history  dependent  phenomena.  However,  closed  form  analytical 
solutions  are  available  for  the  creep  compliance  function  and  the 
recovery  response  if  only  kinematic  hardening  is  considered.  When 
both  types  of  hardening  are  included,  the  differential  model 
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follows  creep  data  on  ice  quite  well,  specifically  those  of  Jacka 
(1984).  Predictions  of  the  ratio  of  transient  (delayed  elastic) 
strain  to  total  strain  agree  qualitatively  with  Sinha's  (1979) 
model  if  grain  size  effects  are  taken  into  account. 

The  multiaxial  generalization  of  the  differential  model 
follows  from  conventional  elasticity  theory  and  from  the  rate 
theory  of  flow.  This  eliminates  the  need  for  an  integral 
formulation  under  variable  loading  histories  or  multiaxial 
loading  and  for  generalizing  the  superposition  assumption  for 
nonlinearly  viscoelastic  materials.  Equations  are  derived  for  an 
orthotropic  model  of  incompressible  flow  and  for  estimating  model 
parameters  from  uniaxial  experimental  data. 

2.  UNIAXIAL  DIFFERENTIAL  MODEL 

Physical  Basis  of  Deformation  Model. —  There  is  general 
agreement,  based  on  theoretical  and  experimental  work,  that  at 
least  two  thermally  activated  deformation  systems,  a  soft  system 
and  a  hard  system,  are  present  during  the  flow  of  fresh-water 
polycrystalline  ice  (Sinha,  1979,  Ashby  and  Duval,  1985.)  They 
may  be  either  grain  boundary  sliding  (with  diffusional 
accommodation)  and  basal  slip  or  basal  slip  and  slip  on  a 
non-basal  plane.  ^  combination  of  these  processes  could  be 
present  as  well. 

Initially,  the  solid  resists  the  applied  stresses  in  an 
elastic  manner  and  then  flow  begins  on  the  soft  and  hard  systems. 
However,  flow,  particularly  on  the  easy  soft  system,  causes  the 
build-up  of  internal  elastic  stresses.  This  may  occur  as  a  result 
of  grain  boundary  sliding  next  to  grains  poorly  aligned  for 
deformation  or  dislocation  pile-ups  at  the  boundaries  of  such 
grains.  Dislocation  pile-ups  at  grain  boundaries  have  been 
observed  in  ice  through  scanning  electron  microscopy  (Sinha, 
1987.)  The  internal  elastic  stresses,  termed  back  or  rest 
stresses ,  resist  flow.  In  addition,  internal  drag  stresses  which 
resist  dislocation  fluxes  are  generated  in  annealed  materials 
undergoing  flow.  The  increase  in  drag  stresses  are  the  outcome  of 
creep  resistant  substructures,  i.e.,  subgrains  and  cells,  formed 
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by  grain  boundary  sliding  or  dislocation  movement  and  of 
dilocation  entanglement,  dipole  formation  and  kink  band  formation 
during  slip  (particularly  on  the  basal  plane.) 

A  detailed  understanding  of  evolving  structural  and  stress 
states  on  the  deformation  of  polycrystalline  materials  is 
unavailable  at  the  present  time.  For  example,  only  recently  has 
an  attempt  been  made  to  model  the  primary  creep  process  resulting 
from  sub-cell  formation  using  sub-cell  size  and  misor ientation  as 
state  variables  (Derby  and  Ashby,  1987.)  However,  it  is  well 
known  that  an  increasing  drag  stress  contributes  to  isotropic 
hardening,  while  an  increasing  rest  stress  contributes  to 
kinematic  hardening.  In  isotropic  hardening,  material  properties 
are  independent  of  the  direction  of  straining.  On  the  other  hand, 
kinematic  hardening  induces  directionally  dependent  material 
properties,  referred  to  as  deformation  or  stress-induced 
anisotropy.  The  Baushinger  effect  in  metals  is  an  example  of 
kinematic  hardening. 

In  this  paper,  the  deformations  resulting  from  the 
interactions  between  the  soft  and  hard  systems  are  decomposed 
into  two  components;  a  transient  flow  component  and  a  steady  flow 
component.  Steady  state  flow,  representing  a  balance  between 
work-hardening  and  recovery,  is  associated  with  viscous 
(irrecoverable)  strains,  isotropic  and  kinematic  hardening 
phenomena  are  active  during  transient  flow  and  give  rise  to 
elastic  strains.  These  strains  are  recoverable  on  unloading  since 
equilibrium  requires  the  internal  elastic  back  stress  to  reduce 
to  zero.  The  time-dependent  elastic  strains  defining  transient 
deformation  represent  the  phenomenon  of  delayed  elasticity  or 
anelasticity . 

Mathematical  Formulation. —  The  governing  equation  for  the 
model  under  uniaxial  conditions  is  obtained  by  expressing  the 
total  strain  rate  as  a  sum  of  its  components,  i.e.: 


c  ■ 


%  ^  "t  s 


(1) 


where  the  three  terms  on  the  right  hand  side,  representing 


5 


instantaneous  elasticity,  transient  flow  and  steady  state  or 
viscous  flow  are  described  in  what  follows. 

The  instantaneous  elastic  strain,  c^ ,  is  related  to  the 
stress,  o,  throuch  the  Young's  modulus,  E,  of  polycrystalline 
ice;  this  relationship  may  be  expressed  in  rate  form  as: 


•  • 

Eg  -  O/E  (2) 

Several  investigators  (see,  e.g..  Gold,  1977)  have  shown  using 
high-frequency  sonic  methods  that  the  Young's  modulus  of 
polycrystalline  fresh-water  ice  varies  in  the  range  of  9-11  GPa, 
with  negligible  temperature  dependence  between  -5®C  and  -45°C. 

The  viscous  strain,  c^,  which  is  associated  with  secondary 
creep  or  steady  flow  conditions,  follows  the  well  Itnown  Norton 
type  power  law  of  Glen  (1955),  i.e.. 


a 


(3) 


where  N  is  the  power  law  index  and  V  is  a  temperature  dependent 
constant  characterized  by  an  Arrhenius  activation  energy  law: 


V  -  exp(Q/NRT) 


(4) 


T  is  the  temperature  in  Kelvin,  is  a  temperature  independent 
constant,  Q  is  the  activation  energy,  and  R  is  the  universal  gas 
constant  equal  to  8.32  J  mol~^  k”^.  The  activation  energy  for 
steady  flow  of  columnar-grained  polycrystalline  ice  has  been 
experimentally  determined  by  Gold  (1973)  to  be  65  KJ  mol”^  for 
temperatures  in  the  range  of  -5®C  to  -40®C.  While  the  activation 
energy  for  pure  single  crystals  does  not  change  with  temperature 
up  to  -0.2*C,  Gold  (1983)  suggests  that  Q  varies  at  the  higher 
temperatures  for  polycrystalline  ice  and  that  at  temperatures 
greater  than  -5®C  it  is  probably  closer  to  100  KJ  mol”^.  Similar 
trends  have  been  observed  by  Barnes  et  al.  (1971.) 


The  transient  strain  rate  is  taken  to  follow  a  Norton 
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where  the  variable  R  represents  the  back  stress  and  B  is  a 
non-dimensional  drag  stress.  Implicit  in  the  formulation  of  Eq. 
(5)  are  the  assumptions  that:  (1)  the  exponent  N  is  the  same  as 
that  for  steady  flow  in  Eq.  (2),  and  (ii)  the  temperature 
dependence  of  the  transient  deformation-rate,  represented  by  the 
parameter  V,  is  given  by  an  Arrhenius  law  with  an  activation 
energy  equal  to  that  for  steady  flow.  For  columnar-grained 
polycrystalline  (fresh-water)  ice  the  former  assumption  can  be 
deduced  from  the  numerical  values  for  parameters  in  Sinha's 
(1978)  time-hardening  model,  and  for  dense  isotropic  polycrystals 
from  the  strain-hardening  model  of  Ashby  and  Duval  (1985).  Sinha 
(1978)  has  also  shown  that  the  activation  energy  for 
transient  flow  is  equal  to  67  KJ  mol~^,  which  agrees  well  with 
Gold's  (1973)  data  for  steady  flow  in  the  same  ^ype  of  ice. 

Evolution  equations  must  be  specified  for  R  and  B  which  are 
history-dependent  variables  representing  transient  flow.  Since 
the  transient  strains  are  elastic  in  nature,  the  time  rates  of 
change  of  the  back  and  drag  stresses  are  linearly  proportional  to 
the  transient  strain  rate.  The  following  equations  are  postulated 
to  descibe  the  evolution  of  R  and  B: 


R 

B 


AE 


sgn 


(6) 

(7) 


The  initial  value  of  R  is  zero  for  an  annealed  material  or  for  a 
material  that  has  recovered  from  prior  loading.  On  the  other  hand 
the  initial  value  of  B,  i.e.,  B^,  may  represent  the  annealed 
state  of  the  material  or  some  level  of  initial  hardening 
introduced  by  pre-straining.  Both  A  and  H  are  temperature 
independent  and  dimensionless  variables. 

Under  creep  loading  R  will  asymptotically  increase  to  a 
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value  equal  to  the  applied  stress,  at  which  point  transient  flow 
will  cease.  In  the  case  of  constant  strain  rate  loading,  R 
approaches  the  steady  state  stress  asymptotically.  The  maximum 
value  of  transient  strain  in  both  these  cases  is  given  by 
c.  _.„“o/AE  when  R  is  zero  initially.  A  value  of  A  less  than  one 

u  f  liio  X 

suggests  that  this  magnitude  is  greater  than  the  instantaneous 
elastic  strain.  For  the  same  loading  conditions,  the  drag  stress 

reaches  a  maximum  value  equal  to  B  +Hc.  _  .  This  constraint  on 

the  maximum  value  of  B  states  that  the  isotropic  resistance  to 
transient  flow  is  not  unbounded;  if  it  is  unbounded  and 
approaches  infinity,  the  material  will  lose  its  ability  to 
undergo  further  flow. 

Under  reversed  or  cyclic  loading  R  will  reverse  or  switch 
back  and  forth  between  positive  and  negative  values,  i.e.,  the 
physical  processes  associated  with  kinematic  hardening  can 
locally  relax  or  move  back  and  forth,  thus  preventing  a  continual 
build-up  which  would  lead  to  considerable  hardening.  The  signum 
function  is  used  in  Eq.  (7)  to  ensure  that  B  has  the  same  effect 
on  material  behavior  under  both  compressive  and  tensile  loadings. 
For  instance,  it  can  be  inferred  from  Eq.  (7)  that  B  >  0  during 
both  tensile  and  compressive  creep  tests,  while  it  is  negative 
during  unloading  in  both  types  of  tests.  The  decrease  in  drag 
stress  during  unloading  indicates  a  decreasing  resistance  to 
grain  boundary  sliding  and  dislocation  fluxes.  This  may  arise 
from  a  spatial  bias  in  the  distribution  of  defects  generated  by 
isotropic  hardening  which  favors  regions  of  high  back  stress 
concentration . 

Equations  (l)-(7)  define  the  governing  differential 
equations  for  the  uniaxial  model.  For  creep  loading,  the 
solutions  of  Eqs.  (l)-(4)  are  trivial.  However,  Eqs.  (5)-(7)  are 
coupled  and  numerical  integration  is  necessary  to  compute  the 
transient  strains  if  both  isotropic  and  kinematic  hardening  are 
present.  If  isotropic  hardening  is  absent,  i.e.,  B  is  a  constant, 
analytical  solutions  can  be  obtained  as  shown  in  a  subsequent 
section.  For  a  general  or  variable  loading  history,  the  governing 
equations  are  all  coupled  and  numerical  integration  is  required. 


Model  Formulation  in  Dimensionless  Variables. —  For  the 
special  cases  of  constant  stress  and  constant  strain  rate 
loading,  Ashby  and  Duval  (1985)  have  suggested  that  unique 
relationships  exist  between  certain  dimensionless  variables.  Such 
relationships  are  predicted  by  the  proposed  model  as  shown  below. 

For  creep  of  polycrystalline  ice  at  constant  applied  stress, 
Ashby  and  Duval  (1985)  have  considered  the  following 
dimensionless  variables  for  strain,  strain  rate,  time,  and  the 


back  stress: 

e  ■  cE/ff  (8) 

t  -  t/t^  (9) 

t  -  tfc^E/o  (10) 

R  -  R/o  ( 11  ) 


Substituting  Eqs.  (8)-(ll)  in  Eqs.  (2),  (3)  and  (5)  yields: 


-  1 


(12) 


S  t  (13) 

and 

fey  -  1  (14) 

R  +  Bfe^.^/*^  -  1  (15) 

In  order  that  Eq.  (7)  also  reduces  to  a  dimensionless  form,  the 
hardening  parameter  H  is  defined  as  HE/c.  The  dimensionless 
evolution  equations  can  then  be  expressed  as: 


R 


(16) 


dt 


sgn 


■dicj.' 
.d  t  . 


(17) 


In  the  above  equations  the  differentiation  is  with  respect  to 
dimensionless  time.  Equations  (12)-(17)  show  that  the  model 
predicts  a  unique  relationship  between  the  dimensionless 
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variables  and  is  independent  of  applied  stress  level  and 
temperature . 

Under  constant  strain  rate  loading  the  model  predicts  that  a 
unique  relationship  exists  between  dimensionless  stress  a  and 
dimensionless  time  t,  independent  of  the  applied  strain  rate  fe, 

O 

and  temperature.  Consider  the  following  dimensionless  variables 
for  stresses,  time,  and  strains,  as  suggested  by  Ashby  and  Duval 
(1985)  : 


a 


t 


c 


e 


min 


min 


R 

R  -  - 

e 

min 


o 


(18) 

(19) 

(20) 


where  is  the  stress  corresponding  to  the  minimum  creep  rate 

given  by  Glen's  power  law,  Eq.  (3); 


a 

min 


V  c 


1/N 


(21) 


Substituting  Eqs.  (18)-(20)  in  Eqs. 


R  -  A  —  [a  e.] 
d?  ^ 


^  -  H 

-  [O  e  l 

sgn 

dt 

dt 

Id  t 

(2),  (3),  (5)-(7)  yields: 

(22) 

(23) 

(24) 

(25) 

(26) 


where  c  indicates  da/dt,  and  similarly  for  and  R. 

Upon  substituting  Eqs . ( 22 )-{ 24 )  in  Eq.  (1)  expressed  in 
dimensionless  form,  the  following  equation  is  obtained: 


N 

a  -  1  - 

0  -  R 

+ 

0  T 

B 

Eqs.  (24)-(27)  can  be  integrated  with  the  initial  conditions  of 
zero  dimensionless  stress  and  transient  strain.  As  steady  state 
is  reached,  i.e.,  the  dimensionless  stress  rate  and  transient 
strain  rate  decay  to  z^'ro,  the  above  equations  show  that  the 
dimensionless  stress  and  transient  strain  tend  to  one  and  1/A, 
respectively.  The  stress  at  steady  state  will  therefore  attain 
the  value  of  o  .  given  by  Eq.  (21).  A  single  master  curve  can  be 
used  to  relate  the  dimensionless  stress  and  time  since  the 
temperature  dependent  constant  V  and  the  applied  strain  rate  have 
been  eliminated  from  the  equations.  Experimental  data  is 
currently  unavailable  for  verifying  the  dimensionless 
relationships  under  constant  strain  rate  loading. 

Closed-Form  Analytical  Solutions  for  Creep  and  Recovery 
Response . —  As  previously  stated,  closed  form  analytical 
solutions  exist  for  creep  and  recovery  response  when  isotropic 
hardening  is  absent,  i.e.,  B  is  a  constant.  These  solutions  are 
valuable  since  they  provide  insights  regarding  the  behavior  of 
the  model.  The  analytical  solutions  are  derived  below. 

In  an  ideal  creep  test  the  stress,  o,  is  applied 
instantaneously  and  the  stress  rate  history  is  a  Dirac  delta 
function,  i{t) ,  with  amplitude  e.  This  history  is  zero  for  all  t 
except  at  t»0  where  it  is  infinity  such  that: 

pt* 

c  S ( t )  «  e  ( 28 ) 

}-• 

for  t*  >  0  and  zero  otherwise.  Consequently,  the  initial  strain 
rate  predicted  by  the  model  is  also  a  Dirac  delta  function,  i.e.. 
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it  is  equal  to  infinity.  The  amplitude  of  this  function  is  o/E, 
which  when  integrated  in  a  manner  similar  to  Eq.  (28)  corresponds 
to  the  instantaneous  elastic  strain.  For  time  incrementally 
greater  than  zero,  the  strain  rate  is  finite  and  equals: 

-  (a/BV)^  +  (a/V)^  (29) 

Equation  (29)  recognizes  that  the  elastic  back  stress  in  Eq.  (5) 
is  equal  to  zero  initially.  Since  the  first  term  of  the  equation 
which  represents  transient  flow  dominates  the  initial  creep 
response,  the  constant  B  will  generally  be  less  than  one. 

The  dimensionless  creep  compliance  function  for  the  model, 

J,  is  the  sum  of  the  dimensionless  elastic,  transient  and  viscous 
strains,  respectively,  i.e. 

j  -  Eg  +  (30) 

The  dimensionless  elastic  and  viscous  strains  are  given  in  Eqs. 
(12)  and  (13).  The  dimensionless  transient  strain  can  be 
analytically  derived  from  Eqs.  (15)  and  (16)  with  a  substitution 
of  variables  approach.  In  particular,  define  a  variable  q  as 
follows : 


q  -  1  -  A  E^.  (31) 

Then, 

q  -  “A  e  ^  (32) 

Substitution  of  Eqs.  (31)  and  (32)  into  Eq.  (15)  and  a  separation 
of  variables  yields: 


(33) 


Integrating  Eq.  (33),  applying  the  initial  condition  of  e^  -  0, 
i.e.,  q  ■  1,  and  substituting  for  q  results  in: 
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-  1/A  -  [  +  A^/B^{N-1 ) t (34) 

Equation  (30)  together  with  Eqs.  (12),  (13)  and  (34)  provide  a 
closed  form  analytical  solution  for  the  dimensionless  creep 
compliance  function.  Also,  by  substituting  Eq.  (34)  in  Eq.  (15), 
the  dimensionless  transient  strain  rate  can  be  expressed  in  terms 
of  dimensionless  time  as: 


[B^  A/B(N-1 


(35) 


Equations  (34)  and  (35)  show  that  the  dimensionless  transient 

N 

strain  and  strain  rate  tend  to  1/A  and  1/B  as  dimensionless  time 


tends  to  infinity  and  zero,  respectively. 

If  creep  recovery  is  allowed  to  occur  at  time  t-t^,  the 
elastic  component  of  the  strain  is  recovered  instantaneously 


while  the  viscous  component  is  irrecoverable  and  remains 


unchanged  with  time.  However,  the  transient  strain  will  decay 
with  time  according  to  the  following  closed  form  analytical 
solution  that  can  be  derived  from  Eqs.  (15)  and  (16)  in  a  manner 
similar  to  Eq .  ( 34  ) : 


(A/B)' 


'(N-l)(t-t^)  J 


1/(1-N) 


(36) 


where  is  the  dimensionless  transient  strain  at  the  time  of 

unloading.  Equation  (36)  shows  that  the  dimensionless  transient 
strain  decays  to  zero  with  dimensionless  time  after  unloading. 


3.  EXPERIMENTAL  VALIDATION  OF  UNIAXIAL  MODEL 

This  section  first  identifies  the  uniaxial  model  parameters 
and  discusses  methods  for  determining  them.  Then,  model 
predictions  under  constant  stress  loading  are  verified  against 
the  experimental  data  of  Jacka  (1984)  and  Sinha  (1978).  The  model 
is  also  compared  with  Sinha's  (1979)  predictions  for  the  relative 
contribution  of  transient  strain  to  the  total  strain'  during 
creep.  To  further  demonstrate  the  capability  of  the  model,  the 
predicted  strain  response  under  a  monotonically  increasing  stress 
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history  is  verified  against  Sinha's  (1981)  test  data.  Finally, 
the  creep  and  recovery  response  of  randomly  oriented  snow  ice  is 
studied  using  Brill  and  Camp's  data  talcen  from  Sinha  (1979). 

Parameter  Identification  and  Estimation. —  The  uniaxial 
model  contains  a  total  of  six  parameters;  E,  N,  V^,  A,  H  and  B^. 
For  single  ice  crystals  and  transversely  isotropic  ice,  five 
independent  elastic  moduli  are  needed  to  describe  elastic 
behavior.  Values  for  these  elastic  moduli  are  available  for 
single  crystals  (see  for  example.  Green  and  MacKinnon,  1956).  The 
value  of  E  for  polycrystalline  isotropic  ice  can  be  estimated 
fairly  well  from  the  elastic  moduli  of  single  crystals  (Gammon  et 
al.,  1983).  Typical  values  of  E  for  isotropic  polycrystalline  ice 
are  given  in  Section  2  of  this  paper. 

Based  on  the  results  of  tests  by  a  number  of  researchers 
carried  out  at  -10“C  in  the  stress  range  0.1  to  2  MPa,  Ashby  and 
Duval  (1985)  have  estimated  the  value  of  N  to  be  three  for  the 
creep  of  isotropic  polycrystals  and  two  for  the  basal  glide  of 
monocrystals.  The  use  of  N-3  for  isotropic  polycrystalline  ice  at 
moderate  stresses  is  supported  by  theoretical  models  which  assume 
dislocation  mobility  as  the  rate-controlling  process  (Baker, 
1982).  Sinha  (1978)  has  also  suggested  the  same  value  for  the 
stress  exponent  in  his  equation  for  the  viscous  creep  of 
polycrystalline  ice. 

The  temperature  independent  constant  and  the  activation 

energy  Q  can  be  estimated  from  creep  data  for  various 

temperatures.  From  the  values  of  the  parameters  used  in  Sinha's 

equation  (1978,  1979),  is  estimated  to  be  6.59  x  10“^  MPa  s^^*^ 

-1  ® 

for  Q*67  KJ  mol 

Under  constant  stress  loading,  the  parameters  A  and  B^ 

determine  the  maximum  value  of  the  transient  strain  (e/AE)  and 

N 

the  initial  transient  strain  rate  (a/B^V)  ,  respectively.  The 
constant  A  can  be  estimated  by  subtracting  the  elastic  strain  and 
the  viscous  strain  from  the  total  strain  when  steady  state  is 
reached.  Since  the  total  recoverable  deformation  is  a/E+e/AE,  the 
fully  relaxed  modulus,  equal  to  the  applied  stress  per  unit 
maximum  recoverable  deformation,  is  given  by  EA/(l-«-A).  This 
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allows  A  to  be  computed  from  a  creep  recovery  test  as  well.  The 
constant  can  be  estimated  from  Eq.  (29).  This  requires 
knowledge  of  the  initial  strain  rate  and  the  constant  V.  The 
latter  can  be  computed  from  Eq.  (4),  but  initial  strain  rates 
derived  from  the  measured  initial  strains  may  be  somewhat 
inaccurate  since  experimental  measurements  of  small  strains  tend 
to  be  unreliable  (Jacka,  1984,  Mellor  and  Cole,  1982).  The 
parameter  H  controls  the  amount  of  isotropic  hardening  at  a  given 
time.  It  can  be  estimated  from  creep  strain  and  strain  rate  data 
using  Eq.  (5).  Having  determined  V,  the  viscous  strain  and  strain 
rate  histories  are  known,  and  the  transient  strains  and  strain 
rates  can  then  be  extracted  from  the  creep  data.  Noting  Eqs.  (6) 
and  (7),  Eq .  (5)  can  be  written  in  the  following  way: 

ln[o  -  AEe^  -  -  In  H  +  In  [Vc^.6^.^/^)  (37) 

The  quantity  on  the  left-hand  side  plotted  against  the  second 
term  on  the  right-hand  side  of  Eq.  (37)  is  a  straight  line  and  H 
can  be  computed  from  its  intercept  with  the  y-axis. 

Comparisons  of  the  model  predictions  with  experimental  data 
in  this  paper  is  based  on  the  following  values  for  N,  E,  and 

Q: 


3 

9500  MPa 
6.59x10“'^  MPa 
67  kJ  mol""^ 


^l/N 


Comparison  of  Model  Predictions  with  Jacka's  Creep  Data. — 
Jacka  (1984)  has  published  results  of  uniaxial  compression  tests 
on  isotropic  polycrystalline  ice  with  a  mean  grain  size  of 
1.7+0. 2  mm.  The  samples  were  tested  under  constant  stress  ranging 
from  0.1  to  1.5  MPa  at  the  following  specific  temperatures:  -5.0, 
-10.6,  -17.8  and  -32.5"C.  Figs.  1,  2  and  3  show  plots  for 
Jacka's  data  (taken  from  Ashby  and  Duval,  1985)  corresponding  to 
c  versus  t,  t  versus  t,  and  t  versus  c,  respectively.  The 
predictions  of  the  model,  obtained  by  solving  Eqs.  (15)-(17)  are 
indicated  by  solid  lines  with  A-0.017,  Bq-0.24  and  H-  0.024. 

Also  shown  are  the  model  predictions  for  no  isotropic  hardening. 
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i.e.,  solutions  provided  by  Eqs.  (12)-(14)  and  Eqs.  (34)-(35). 
Parameters  used  for  generating  these  curves  are  A-0.033  and 
B-0.402.  Note  that  for  experimental  data  plotted  in  dimensionless 
form,  the  m.odel  predictions  using  Eqs.  (12)-(17)  as  well  as  Eqs. 
(34)-(35)  are  independent  of  E,  and  Q.  Ashby  and  Duval  (1985) 
have  modified  Sinha's  equation  for  creep  to  a  form  which 
satisfies  the  dimensional  requirements.  The  predictions  of  the 
modified  equation  are  also  shown  in  the  figures.  In  referring  to 
Jacka's  data  it  is  understood  that  all  variables  are  normalized 
and  the  word  dimensionless  is  dropped  when  referring  to  them. 

The  solid  lines  show  that  agreement  between  model 
predictions  and  data  is  good  when  strain  rate  is  plotted  against 
time  or  strain  (Figs.  1  and  2).  The  predicted  master  curve  in  the 
strain  versus  time  plot  (Fig.  3)  represents  the  data  well,  but 
at  small  times  the  predicted  strains  somewhat  overestimate  the 
experimental  data.  With  no  isotropic  hardening,  the  strain  rates 
are  underestimated  while  the  strains  agree  well  with  data.  On  the 
other  hand,  the  modified  equation  provides  a  good  prediction  of 
the  strain  rate  versus  time  response  (Fig.  1),  but  it  over¬ 
predicts  the  initial  strains  (Figs.  2  and  3)  in  spite  of  a  factor 
of  two  reduction  in  the  value  of  the  parameter  (parameter  A  in 
Ashby  and  Duval's  paper)  which  equals  the  maximum  transient 
strain . 

Comparison  of  Model  Predictions  with  Sinha's  Creep  Data. — 
Sinha  (1978)  has  conducted  tests  on  the  creep  behavior  of 
transversely  isotropic  columnar-grained  ice  (S-2  ice)  with  an 
average  grain  diameter  of  3  mm.  The  tests  were  conducted  in  the 
temperature  range  of  -9.9  to  -41®C  under  a  uniaxial  compressive 
load  of  0.49  MPa  acting  in  the  plane  of  transverse  isotropy. 

Based  on  the  observations  that  the  activation  energy  for  both 
viscous  flow  and  transient  deformation  appears  to  be  equal  and 
that  Young's  modulus  is  relatively  independent  of  temperature, 
Sinha  (1978)  postulated  that  the  time  dependence  of  the  strain  at 
one  temperature  can  be  obtained  by  shifting  the  measured 
dependence  at  another  temperature  along  the  time  scale  using  a 
shift  function.  Figure  4  shows  the  creep  strains  obtained  at 
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various  temperatures  shifted  to  a  reference  temperature  of  -10®C. 

The  solid  line  indicates  the  model  prediction  with  A-0.33, 

B  -0.058  and  H-0.63.  The  value  of  A  is  identical  to  that  obtained 
o 

by  Sinha  ( 1978  )  . 

The  agreement  between  the  experimental  data  and  theoretical 
results  is  very  good.  Notice  also  that  the  values  of  A,  H  and 
used  for  Sinha's  and  Jacka's  data  are  different,  reflecting 
differences  in  the  ice  types  that  were  tested,  i.e.,  isotropic 
and  granular  versus  transversely-isotropic  and  columnar-grained, 
and  the  average  diameters  of  ice  grains.  Such  modifications  to 
parameter  values  are  also  needed  for  Sinha's  equation.  For 
example,  the  parameter  corresponding  to  A  was  determined  to  be 
1/3  from  Sinha's  tests  on  ice  with  a  grain  size  of  3  mm,  but 
values  of  1/70  and  1/35  were  found  to  be  suitable  for  Jacka's 
data  (Ashby  and  Duval,  1985). 

Model  Prediction  of  Ratio  of  Transient  to  Total  Strain. — 

The  parameter  AE  can  be  interpreted  as  an  anelastic  modulus  (not 
to  be  confused  with  relaxed  modulus),  while  BV  represents  the 
resistance  to  transient  flow.  If  transient  deformation  is  related 
to  grain  size  as  postulated  by  Sinha  (1979),  then  the  parameters 
A  and  H  will  depend  on  grain  size.  Sinha's  (1979)  model  considers 
both  transient  strain  and  strain  rate  to  be  inversely 
proportional  to  grain  size.  For  consistency  with  this 
formulation,  it  is  necessary  for  the  model  parameters  to  be 
related  to  the  grain  size  d  as  follows: 


A  -  d/A' 

B  -  B 
0  o 


(38) 

(39) 


H  -  H'  (40) 

where  A',  B^'  and  H'  are  grain  size  independent  material 

parameters.  The  values  of  these  latter  parameters  are  calculated 

from  Eqs.  (38)-(40)  respectively.  For  the  previously  determined 

values  of  A,  B  and  H^(viz.,  0.33,  0.058  and  0.63)  A'-9  mm, 

B^'-0.04  mm"^-^^,  and  H'-0.3  mm”^/^. 
o 

The  analytical  solutions  for  the  case  of  no  isotropic 
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hardening  are  useful  for  inferring  some  of  the  characteristics  of 
the  complete  model.  It  is  apparent  from  Eqs.  (14)  and  (35)  that 
the  ratio  of  the  (dimensionless)  transient  strain  rate  to  the 
viscous  strain  rate  decreases  with  increase  in  grain  size.  The 
ratio  also  decreases  with  dimensionless  time  (Eq.  35).  According 
to  the  definitions  of  dimensionless  time  (Eq.  (10))  and  viscous 
strain  rate  (Eq.  (3)),  the  ratio  must  also  decrease  with  increase 
in  applied  stress.  These  trends  are  in  agreement  with  predictions 
of  Sinha's  equation. 

The  predictions  of  the  proposed  model  and  Sinha's  equation 
with  regard  to  the  relative  contribution  of  transient  strain  to 
total  strain  are  compared  below.  Let  the  ratio  of  the  transient 
strain  to  the  total  strain,  r,  be  defined  as  follows: 


c  +  t  +  1 

where  the  dimensionless  variables  have  been  defined  previously. 

By  solving  Eqs.  (15)-(17)  the  strain  dependence  of  y  under  a 
constant  stress  load  of  1  MPa  at  -10*C  for  various  grain  sizes  is 
predicted  as  shown  in  Fig.  5.  The  important  features  predicted  by 
Sinha's  equation  such  as  the  increasing  value  of  y  with 
decreasing  d,  the  occurrence  of  maximum  y  at  small  strains,  the 
gradual  shift  of  the  maxima  towards  larger  strains  with 
decreasing  d,  the  gradual  decrease  in  y  with  increase  in  strain 
after  the  peak  is  passed,  and  the  decreasing  effect  of  d  on  y  at 
large  strains  are  also  observed  in  Fig.  5,  although  the  numerical 
values  are  different. 

Furthermore,  since  the  relationship  between  the 
dimensionless  transient  strain  and  time  is  independent  of 
temperature  (Eq.  (34)),  it  can  be  deduced  from  Eq.  (41)  that  the 
evolution  of  y  with  dimensionless  strain  for  a  given  grain  size 
is  unique,  i.e.,  independent  of  both  temperature  and  stress 
level.  Recalling  that  dimensionless  strain  is  equal  to  cE/c  and 
if  E  does  not  change  appreciably  with  temperature,  then  for  a 
given  grain  size,  the  evolution  of  y  with  strain  (stress)  itself 


is  independent  of  temperature  but  not  of  stress  (strain).  This  is 
also  predicted  by  Sinha's  equation. 

Prediction  of  Model  Response  Under  Monotonically  Increasing 
Stress . —  The  rate  sensitivity  of  the  compressive  strength  of 
columnar-grained  ice  under  constant  cross-head  displacement  rates 
has  been  investigated  by  Sinha  (1981).  It  was  shown  that  the 
results  are  representative  of  the  constant  stress  rate  rather 
than  the  constant  strain  rate  condition.  A  numerical  integration 
method,  based  on  a  generalized  creep  equation  and  the  principle 
of  superposition,  was  developed  by  Sinha  (1983)  to  predict  the 
evolution  of  strain  corresponding  to  a  given  stress  history. 

For  the  proposed  model,  the  strain  response  can  be  obtained 
by  numerically  integrating  Eqs.  (l)-(7).  In  this  example,  the 
actual  stress-time  history  (not  the  constant  stress  rate 
idealization)  is  taken  as  input  and  is  known  from  Sinha's  (1981) 
tests  on  ice  with  an  average  grain  size  of  4.5  mm.  The  tests  were 
carried  out  at  -10®C  under  a  constant  cross-head  displacement 
rate  of  1.25x10“^  cm/s.  The  values  of  the  hardening  parameters 
are  determined  from  Eqs.  (38)-(40)  for  the  given  grain  size  and 
previously  determined  values  of  the  grain  size  independent 
parameters . 

The  stress-time  data  is  presented  in  the  upper  curve  of  Fig. 
6b,  while  the  lower  curve  shows  the  predicted  strain-time 
response  superimposed  on  the  test  data.  The  agreement  between 
theory  and  experiment  is  quite  good,  given  that  the  parameter 
values  which  were  determined  from  a  different  data  set  are 
unchanged.  Figure  6a  shows  that  when  stress  is  plotted  against 
strain,  a  very  good  representation  of  the  data  is  obtained.  It  is 
possible  to  conclude  from  this  figure  that  the  predictions  of  the 
proposed  model  under  monotonically  increasing  stress  compare  well 
with  experimental  data. 

Comparison  of  Model  Predictions  with  the  Creep  and  Recovery 
Data  of  Brill  and  Camp. —  Figure  (7)  shows  creep  and  recovery 
data  for  tests  conducted  on  randomly  oriented  snow  ice  by  Brill 
and  Camp  (reproduced  in  Sinha,  1979.)  The  three  sets  of  data 
refer  to  tests  carried  out  under  the  following  conditions:  curve 
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(a)  at  -5®C  and  0.232  MPa,  curve  (b)  at  -5®C  and  0.125  MPa,  and 
curve  (c)  at  -10®C  and  0.238  MPa.  The  model  predictions,  shown  in 
solid  lines,  are  generated  with  A'«6.5  mm,  B^'-O.ll  mm  '  and 

H'-O.Ol  mm”^^^.  The  grain  sizes  used  for  curves  (a),  (b)  and  (c) 
are  2  mm,  2.3  mm,  and  1.5  mm  respectively,  which  are  almost 
identical  to  the  values  determined  by  Sinha  (1979).  Differences 
in  the  hardening  parameters  reflect  the  difference  in  ice  types, 
i.e.,  transversely  isotropic  and  columnar-grained  versus 
isotropic  and  granular  snow  ice.  The  agreement  between  model 
predictions  and  test  data  is  quite  good,  given  that  the 
measurement  of  strain  recovery  in  ice  shows  large  scatter  (Sinha, 
1982.  ) 

A  major  difference  exists  between  Sinha's  recovery  model  and 
the  present  formulation.  The  former  can  result  in  a  decrease  of 
the  permanent/irrecoverable  viscous  strain  and  eventually  lead  to 
reversed  strain  (e.g.,  an  elongation  or  tensile  strain  due  to 
recovery  from  compressive  creep).  This  is  due  to  the  particular 
form  of  the  superposition  principle  adopted,  in  which  the  elastic 
and  the  transient  strains  resulting  from  the  stress  drop  are 
subtracted  from  the  total  strain  at  unloading.  Recovery  is  thus 
the  mirror  image  of  the  transient  term  in  the  equation  for 
loading  and  as  time  increases  it  can  exceed  the  transient  creep 
strain  at  the  instant  of  unloading.  In  order  to  overcome  the 
problem,  the  superposition  principle  is  not  used  when  the 
predicted  strain  during  recovery  becomes  less  than  the 
accumulated  viscous  strain  and  the  strain  is  kept  fixed 
thereafter  at  a  value  equal  to  that  of  the  viscous  strain  at 
unloading.  The  proposed  theory  does  not  suffer  from  this  modeling 
limitation  since  the  values  of  R  and  6  decrease  during  unloading, 
reflecting  creep  recovery. 

4.  MULTIAXIAL  MODEL  FORMULATION 

Natural  ice  has  very  complex  crystalline  and  stratigraphic 
structures,  and  generally  cannot  be  considered  as  an  isotropic 
material.  For  example,  columnar  fresh-water  ice  may  have  two 
sources  of  anisotropy:  (a)  the  c-axis  may  be  oriented 
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perpendicular  to  the  axis  of  crystal  growth,  and  (b)  the  c-axes 
of  different  crystals  may  show  preferred  orientation  in  the  plane 
on  which  they  lie.  According  to  the  classification  of  Michel  and 
Ramsc^er  (1971),  the  first  source  of  anisotropy  is  exhibited  by 
S2  ice  while  both  types  of  anisotropy  are  present  in  S3  ice. 

The  anisotropy  of  ice  strongly  influences  its  mechanical 
behavior.  Carter  and  Michel  (1971)  have  tested  S2  ice  at  -10®C 
under  constant  strain  rate  loading  conditions.  They  find  that  the 
first  source  of  anisotropy  leads  to  a  vertical  to  horizontal 
maximum  stress  or  strength  ratio  of  about  two.  Information  on  the 
effect  of  the  second  source  of  anisotropy  on  the  strength  ratio 
of  freshwater  ice  is  currently  unavailable,  although  data  for  sea 
ice  indicates  the  following  strength  ratios:  (a)  0.25-0.60  for 
strength  at  a  45  degree  azimuthal  angle  to  that  along  the  c-axes, 
and  (b)  0.50-0.95  for  strength  at  a  90  degree  angle  to  that  along 
the  c-axis  (Peyton,  1968;  Vittoratos,  1979;  Wang,  1979; 
Richter-Menge  et  al . ,  1985). 

Theoretical  formulations  which  account  for  anisotropy  in  ice 
with  a  transversely  isotropic  model  have  been  developed  by 
Reinicke  and  Ralston  (1977)  and  by  Vivatrat  and  Chen  (1985).  The 
former  model  is  based  on  plasticity  theory  and  considers  ice  to 
be  a  pressure  sensitive  material  as  well.  The  latter  is  a 
pressure  insensitive,  elastic  -  power  law  creep  formulation. 

The  development  presented  here  is  based  on  an  orthotropic 
generalization  (i.e.,  the  general  case  of  a  material  having  three 
orthogonal  planes  of  symmetry)  of  the  proposed  uniaxial  model 
which  accounts  for  both  transient  and  steady  state  flow  in  ice. 
The  transversely  isotropic  and  isotropic  formulations  are  special 
cases  of  the  orthotropic  generalization. 

Conceptual  Framewor):  and  Constraint  Conditions. —  The 
three-dimensional  generalization  of  the  model  follows  naturally 
from  the  uniaxial  formulation,  i.e.,  it  is  based  on  strain 
decomposition,  linear  elasticity  and  the  rate  theory  of  flow. 
Constitutive  relations  are  derived  for  each  mechanism  of 
deformation  in  the  model,  resulting  in  the  orthotropic  equivalent 
of  Eqs .  ( 1 )-( 7  ) . 
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To  model  orthotropic  elasticity,  the  classical  formulation 
from  elasticity  theory  is  adopted  (see,  for  example,  the  book  by 
Lekhnitskii,  1963).  Compressibility  of  ice  deformation  is 
implicitly  contained  in  linear  elasticity  where  the  Poisson's 
ratios  are  less  than  0.5,  while  the  transient  and  viscous 
deformation-rate  mechanisms  are  assumed  to  be  incompressible 
(Palmer,  1967,  Sinha,  1987'.  The  orthotropic  generalization  of 
the  viscous  deformation-rate  mechanism  is  derived  from  the  rate 
theory  of  flow  by  applying  the  normality  principle  to  a  scalar 
valued  flow  potential  expressed  in  terms  of  an  equivalent  stress 
measure  for  incompressible  orthotropic  materials.  Similarly,  the 
derivation  of  the  orthotropic  constitutive  relations  for 
transient  flow  are  based  on  the  normality  of  the  stress 
difference  vector  ^  -  R  to  a  scalar  valued  flow  potential 
expressed  in  terms  of  an  equivalent  stress  difference  measure. 
Evolution  equations  for  the  back  stress  vector  R  and  the 
equivalent  drag  stress  follow  from  the  uniaxial  equations. 

The  total,  elastic,  viscous  and  transient  strain  rate 
vectors  must  obey  the  following  constraint  condition: 

•  «  •  • 

c  -  Cg  +  Cv  +  It  (42) 

where  the  strain  rates  are  in  engineering  notation,  for 
example,  fe  -  ( ^he  superscript  T 
denotes  the  transpose  of  vectors  and  matrices.  For  convenience, 
the  stress  difference  c-R  is  denoted  by  the  symbol  i.e.: 

£(j  “  £  ”  B  (43) 

The  stress  vectors  are  also  expressed  in  engineering  notation. 

Orthotropic  Model  of  Elasticity. —  The  constitutive  relation 
between  the  elastic  strain  and  the  stress  is  described  in  rate 
form  as: 

-  C  ff  (44) 

where  C  is  the  compliance  matrix  for  a  linearly  elastic  but 


22 


orthotropic  material.  Values  of  the  Young's  moduli,  Poisson's 
ratios  and  shear  moduli  for  orthotropic  and  transversely 
isotropic  polycrystalline  ice  are  not  readily  available.  However, 
engineering  approximations  involving  a  weighted  average  of  the 
five  elastic  constants  for  single  crystals  have  been  developed 
(Gammon  et  al.,  1983,  Ashton,  1986).  The  Poisson's  ratio  for 
isotropic  polcrystalline  ice  is  approximately  0.3  (Gold,  1977). 

Orthotropic  Model  of  Viscous  Flow. —  To  derive  the 
relationship  between  the  viscous  strain  rate  vector  and  the 
stress  vector  a,  an  equivalent  stress  measure  generalized  for 
pressure  insensitive  orthotropic  materials,  i.e.,  with  identical 
behavior  in  tension  and  compression,  is  defined: 


I  ^2  ®2,  ,2  ®3,  ,2 

eq  6  [3  W  3  VY  zz  3  zz  xx ' 


(45) 


with  e  chosen  to  be  (a.+a.,)  so  that  when  the  stress 

i  z  e  yy 

components  are  described  by  the  vector  a  ■  (0  0000)  , 

i.e.,  the  y-axis  is  chosen  as  the  reference  direction.  Equation 
(45)  is  similar  in  form  to  that  used  by  Hill  (1950)  for  metal 
plasticity  and  may  be  expressed  in  compact  form  using  matrix 
notation  as: 


3/6  G  a 


where 


G 


3  3  3 

^1*^2  "®2 
3  3 

^2-^°3 
3 

SYMMETRIC  ^®4 


0 


(46) 


(47) 
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The  viscous  strainrate  vector  can  now  be  related  to  the  stress 
vector  by  defining  a  scalar  valued  viscous  flow  potential 
function ; 

N+1 

■  a  (  48  ) 

^  N+1 

which  obeys  the  normality  principle: 


(49) 


The  parameter  a  in  Eq.  (48)  is  a  constant  associated  with  the 

power  law  for  uniaxial  loading  in  the  y-direction;  it  is 

N 

equivalent  to  the  quantity  1/V  in  Eq.  (3).  Combining  Eqs. 
(46)-(49)  yields  the  desired  relationship: 


fy  -  X  S*  (50) 

where 

X  -  3/e  a  (51) 

and 

S*  -  G  ff  (52) 


Note  that  S*  is  a  pseudo-deviatoric  stress  vector  for  orthotropic 

materials.  If  a^^  to  a^  ■  1,  S*  reduces  to  the  conventional 

deviatoric  stress  vector,  reduces  to  the  conventional 

eq 

equivalent  stress  measure  for  isotropic  materials,  and  Eq.  (50) 
becomes  the  well  known  three-dimensional  generalization  of  the 
power  law  for  creep  of  isotropic  materials,  as  presented  by 
Palmer  (1967)  for  glacier  flow. 

Using  the  hypothesis  of  energy  equivalence,  the  relationship 
between  the  equivalent  stress  defined  in  Eq.  (45)  and  an 
equivalent  strain  rate  measure  can  be  established.  The  rate  of 
dissipation  of  energy  per  unit  volume,  P,  is  given  by: 


(53) 
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Application  of  this  hypothesis  yields: 


a  6  “  o  t 

-  -V  eq  v,eq 


(54) 


where  c..  is  the  equivalent  viscous  strainrate.  The  viscous 
V ,  eq 

strain  rate  vector  in  Eq.  (54)  can  be  eliminated  using  Eqs.  (50), 
(52)  and  (46)  in  succession  to  yield: 


v,eq 


(55) 


Given  the  equivalent  stress  measure,  Eq.  (55)  can  be  used  to 
compute  the  equivalent  viscous  strain  rate.  Alternatively,  an 
explicit  expression  can  be  derived  by  first  eliminating  {e 


Oyy)^ 


,  in  Eq.  (45)  through  the  use  of  Eq.  (50)  and  then 


substituting  the  resulting  expression  for  in  Eq .  (55).  The 

final  expression  can  be  expressed  in  compact  notation  as  follows: 


’  2  ’  T  ’ 

e  »  6/3  c  H  c 

V ,  eq  iv  -  ~v 


where  the  transformation  matrix  H  is  given  by: 
p(a^+a3)a2^  -3aia2a3  -3a^a2a3 


(56) 


3(  aj^+a2  )a3' 


-3a  j^a2a3 


3(a2+a3)aj^‘ 


(57) 


SYMMETRIC 


with  a*  «  a2a34a3a2’'’a2a2  •  It  is  apparent  that  Eq.  (56)  can  be 
reduced  to  the  conventional  equivalent  strain  rate  measure  for 
isotropic  materials  if  a^  to  a^  •  1.  Moreover,  when  loading  is  in 
the  reference  direction,  and  Eq.  (55)  reduces  to 

uniaxial  loading  in  the  reference  direction. 

Orthotropic  Model  of  Transient  Flow. —  The  orthotropic 
generalization  of  the  transient  deformation  is  based  on  the 
assumption  of  flow  incompressibility.  Although  this  may  not  be 
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strictly  true,  Sinha  (1987)  has  argued  that  transient  flow  does 
not  change  the  volume  appreciably  and  that  the  assumption  is 
valid.  This  assumption  is  implicitly  made  for  metals  (Hart, 

1976).  A  second  assumption  is  that  the  orthotropy  is  described  by 
the  same  set  of  parameters,  i.e.,  a^  through  a^. 

The  model  accounts  for  isotropic  hardening  as  well  as 
kinematic  or  directional  hardening,  which  leads  to  subsequent 
deformation  or  stress-induced  anisotropy  (as  opposed  to  material 
or  texture  anisotropy).  The  relationship  between  the  transient 
strain  rate  and  the  stress  vector  can  be  derived  from  the 
normality  of  to  a  scalar  valued  transient  flow  potential 
function.  Following  the  procedure  used  to  derive  Eqs.  (50)-(52) 
yields : 


3/6  a/B. 


N-1  -  * 
°d,eq  -d 


(58) 


where 


* 

**  £  £(5  “  £  £  ~  £  5 


(59) 


is  the  equivalent  nondimensional  drag  stress,  and  a, 
eq  ^  ^  dfeq 

3/6  £jj  G  £^.  To  complete  the  multiaxial  formulation,  the 

evolution  equations  for  R  and  B  ^  as  well  as  the  value  cf  the 

eqivalent  stress  difference  measure  are  required. 

For  consistency  with  the  incompressiblilty  constraint  on 

transient  flow  and  the  elastic  nature  of  back  stresses,  it  is 

necessary  to  define  a  scalar  valued  flow  potential  in  terms  of 

2  T 

the  equivalent  back  stress,  »3/6  R  £ 

♦s  - 

where  b  equals  1/AE  in  the  reference  direction.  The  transient 
strain  can  then  be  related  to  the  pseudo-deviatoric  back  stress 
vector  by  imposing  normality: 


(61) 
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where 

Ir*  -  5  ^ 


(62) 


The  evolution  equation  for  the  bade  stress  vector  is  the  time 
derivative  of  Eq.  (61),  i.e.. 


Sr 


-  G  R  -  e/(3b) 


(63) 


The  equivalent  non-dimensional  drag  stress  is  given  by: 


®eq  “ 


■t,eq 


sgn 


dc 


t  ,eq 


dt 


(64) 


where  c  equals  H  in  the  reference  direction.  Both  the  equivalent 
transient  strain  rate  and  strain  measures  in  Eq.  (64)  can  be 
obtained  using  the  transformation  matrix  H,  derived  in  Eq.  (57) 
for  the  equivalent  viscous  strain  rate. 

The  equivalent  stress  difference  can  be  expressed  as  a 

2 

function  of  the  equivalent  transient  strain.  Noting  that  o. 

if  T  *  ^  f T 

can  be  defined  in  terms  of  S,  as  3/|5  S.  ,  and  substituting 

S  -Sp  for  and  o-R  for  yields: 


‘'d,eq^  “  ~  ^  Ir* 

*  ★ 

£  and,  consequently,  S  may  be  considered  as  given,  while  Sp  may 

be  computed  by  integrating  Eq.  (63).  Substitution  of  Eq.  (61)  in 

T  T 

the  last  term  of  Eq,  (65)  yields  (e/3b)  R  where  R  equals 
twice  the  elastic  strain  energy  stored  in  the  material.  Based  on 
an  equivalence  in  the  rate  of  stored  elastic  strain  energy,  it 
follows  that: 


St  - 


■t,eq 


(66) 


The  equivalent  uniaxial  relationship  between  R  and  c.  is 

eq  t , eq 

given  by: 


1/b  c 


t  ,eq 


(67) 


since  b-l/AE.  Equation  (67)  and  the  right-hand-side  of  Ec.  (66) 

T 

are  integrated  simultaneously  with  respect  to  time  to  yield  R  c 
On  substitution  in  Eq.  (65),  the  following  result  is  obtained: 


"d,eq"  -  3/e 


2/b  /cj  * 


The  second  term  in  this  equation  is  obtained  by  substituting  Eq 
(61)  in  the  corresponding  term  of  Eq.  (65). 

Equations  (42),  (44),  (50),  (58),  (63)  and  (64)  are  the 
orthotropic  counterparts  of  Eqs.  (l)-(3)  and  (5)-(7).  They  form 
the  governing  equations  that  can  be  integrated  numerically  to 
predict  the  model  response  under  variable  loading  histories 
involving  multiaxial  states  of  stress. 

5.  EXPERIMENTAL  VALIDATION  OF  MULTIAXIAL  MODEL 

Estimation  of  Orthotrooic  Model  Parameters. —  The 


orthotropic  model  parameters  can  be  estimated  from  experimental 
data  under  steady  viscous  flow  conditions.  Five  uniaxial 
(compression)  tests  are  required  to  obtain  the  five  parameters  82 
through  a^  since  a^^  can  be  set  to  one  without  loss  of  generality. 
In  a  comprehensive  paper  reviewing  the  constants  used  in  Glen's 
power  law  for  polycrystalline  glacier  ice,  Hooke  (1981)  has 
concluded  that  in  the  absence  of  experimental  evidence  to  the 
contrary,  a  value  of  three  for  the  power  law  index  N  is 
reasonable,  irrespective  of  the  "structural  state",  e.g.,  fabric 
and  grain  size.  The  effect  of  the  structural  state  is  then 
accounted  for  by  changing  the  "viscosity"  parameter  (V  in  the 
present  model).  This  is  the  approach  adopted  here,  in  which  N  is 
three  and  the  initial  texture  or  material  anisotropy  is  accounted 
for  through  the  use  of  an  appropriate  equivalent  stress  measure. 
Under  uniaxial  loading  in  any  specific  direction,  the  viscosity 
parameter  relating  viscous  strain  rate  and  stress  in  the 
specified  direction  is  provided  by  Eq.  (55),  the  definition  for 
the  equivalent  stress  in  Eq.  (45),  and  the  definition  for  the 
equivalent  viscous  strain  rate  in  Eq.  (56). 
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The  x-axis  is  taken  to  be  normal  to  the  ice  sheet  which  is 
defined  by  the  y-z  plane.  The  c-axes  of  the  ice  crystals  are 
assumed  to  lie  in  the  y-z  plane  and  are  aligned  in  the 
y-direction.  The  tests  are  conducted  in  three  orthogonal 
directions  y,  x,  and  z  respectively,  and  along  the  three  45“  axes 
on  the  y-z,  x-y,  and  z-x  planes  respectively.  Let  |3j^  to  Pg 
represent  the  experimentally  determined  ratios  of  the  maximum 
stresses  (strengths)  for  the  latter  five  tests,  respectively,  to 
the  maximum  stress  in  the  reference  y-direction  for  tests 
conducted  at  the  same  constant  strain  rate.  In  the  case  of  creep 
tests,  the  g's  represent  inverse  ratios  of  the  corresponding 
minimum  strain  rates  raised  to  the  power  of  1/N.  The  parameters 
3^2  to  t)e  determined  from  the  following  equations  (see 

Appendix  A  for  derivations): 


^1 

-  ft  ” 
^2 

(1-6^") 

®2 

ft  ^ 
^1 

- 

(1  +  3;^") 

\  vy  f 

(1-33") 

(  1  C\  \ 

^3 

^1 

-  ft  " 

(1  +  3^^^) 

\  f  V  ) 

- 

3/6 

[434“" 

- 

(71) 

“5 

- 

3/6 

(433"" 

- 

(72) 

^6 

■i 

3/6 

(435“" 

-  n 

(73) 

where  n-2N/(N+l).  Typical  values  for  lie  between  2-5.  While 
the  values  of  the  constants  @2  to  3^  are  not  generally  available 
in  the  literature  for  pure  polycrystalline  ice,  they  may  be 
estimated  from  the  sea  ice  data  referred  to  in  the  beginning  of 
Section  4. 

For  a  transversely  isotropic  material,  i.e.,  isotropy  in  the 
y-z  plane,  32“®3“1  P4“®5*  As  a  result,  a^-a^-l,  84=3^,  the 

parameters  a2  and  a^  are  functions  of  only  3^^ ,  while  a^  depends 
on  both  3^  and  3^.  Only  two  uniaxial  tests  are  required  to  obtain 
3j^  and  3^:  one  in  the  x-direction  and  one  along  the  45“  axis  on 
the  x-y  or  z-x  planes. 
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Model  Predictions  Under  Steady  State  Plane-Strain  and 
Triaxial  Compressive  Loadings. —  Experimental  data  on  the  pure 
flow  (both  transient  and  steady  state)  of  polycrystalline  ice 
under  multiaxial  states  of  stress  is  unavailable,  although  the 
incompressible  and  istropic  power  law  of  Palmer  (1967)  is  widely 
used  to  describe  the  deformation  of  glacier  ice  under  such 
stresses.  In  spite  of  data  limitations,  an  attempt  is  made  here 
to  evaluate  model  predictions  under  steady  flow  conditions. 

Frederking  '1977)  has  conducted  plane  strain  uniaxial 

compression  tests  on  columnar-grained  transversely  isotropic 

freshwater  ice.  For  his  type  A  tests,  strains  in  the  z-direction 

are  constrained  to  zero  and  stresses  are  applied  in  the 

y-direction.  At  steady  state  where  the  power  law  orthotropic 

formulation  suffices,  the  ratio  r  of  the  plane  strain  stress  to 

z 

the  unconfined  stress  at  the  same  strain  rate  is  directly  related 
to  by  the  following  equation  (see  Appendix  B) 


2n  -,1/n 


(74) 


"  L46^"  -  IJ 

The  equation  predicts  r  to  vary  between  2. 1-5.1  for 

z 

experimentally  observed  values  of  ranging  from  2  to  5,  and 
N-3.  This  is  consistent  with  Frtderking's  experimental 
observations  which  were  close  to  2  at  high  strainrates  and  to  5 
at  low  strain  rates.  In  his  type  B  tests,  strains  in  the 
x-direction  are  constrained  to  zero  while  stresses  are  again 
applied  in  the  y-direction.  In  this  case,  the  ratio  is  given 

by: 


(75) 


Since  is  generally  greater  than  one,  will  be  less  than 
approximately  1.2  for  N»3.  For  typical  values  of  ,  the 
predicted  values  of  lies  between  1.01  to  1.06.  This  is 
consistent  with  Frederking' s  experiments  which  showed  negligible 
influence  of  x-direction  confinement  on  stresses.  Although  the 
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derivation  in  Appendix  B  determines  and  as  the  ratios  of 
two  steady  state  stresses  resulting  from  viscous  flow, 
considerable  damage  occurred  in  Frederking's  tests.  The  accuracy 
of  the  predictions  are  interesting  nonetheless.  This  probably 
occurs  because  both  the  unconfined  and  partially  confined 
strengths  are  reduced  by  damage,  and  the  resulting  effect  on 
strength  ratios  is  less  significant. 

According  to  the  orthotropic  model,  the  ratio  of  the 
maximum  axial  stress  with  a  confining  pressure  equal  to  t  times 
the  axial  stress  to  the  maximum  axial  stress  in  the  unconfined 
state  at  the  same  strain  rate  should  be  given  by  (see  Appendix 
C)  : 

1 

-  -  (76) 

1-T 

The  shear  stress  (i.e.  axial  stress  minus  radial  stress) 
normalized  by  the  unconfined  stress  is  independent  of  t  or 
confining  pressure  for  the  model  and  equal  to  one.  The  triaxial 
behavior  of  pure  (non-saline)  polycrystalline  ice  has  been 
studied  by  Jones  (1978).  His  tests,  which  were  performed  at 
strain  rates  of  lO"®  to  5x10”^  s-1,  indicate  up  to  a  factor  of 
two  increase  in  shear  stress  due  to  confining  pressure.  Nadreau 
and  Michel  (1986)  have  reported  triaxial  tests  on  freshwater, 
iceberg  and  saline  ice,  and  their  results  confirm  that  shear 
strength  increases  with  confining  pressure  and  strain  rate.  The 
pressure  and  strain  rate  sensitivity  of  damage  in  ice,  which 
causes  the  increase  in  shear  strength  with  confining  pressure  and 
strain  rate,  is  examined  in  the  forthcoming  paper  by  Wu  and  Shyam 
Sunder  ( 1988 ) . 

6.  CONCLUSIONS 

This  paper  presents  a  multiaxial  differential  flow  law  for 
polycrystalline  ice  which  attempts  to  model  the  underlying 
physical  deformation  mechanisms  active  in  the  material. 
Instantaneous  elasticity  is  modeled  by  the  classical  theory  of 
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linear  elasticity,  while  the  steady  viscous  def ormation~rate 
mechanism  is  described  by  Glen's  power  law.  On  the  other  hand, 
the  transient  deformation-rate  mechanism  is  modeled  by  the 
interaction  between  the  soft  and  hard  deformation  systems  which 
gives  rise  to  an  internal  drag  stress  and  a  back  stress. 
Increasing  drag  and  back  stresses  are  associated  with  the 
phenomena  of  isotropic  and  kinematic  hardening,  respectively. 
Dimensional  requirements  identified  by  Ashby  and  Duval  (1985)  are 
satisfied  by  the  model. 

The  multiaxial  generalization  follows  from  conventional 
elasticity  theory  and  from  the  rate  theory  of  flow  for  the 
viscous  and  transient  deformation-rates.  The  rate  theory  assumes 
normality  of  the  deformation-rate  to  a  scalar  valued  flow 
potential  expressed  in  terms  of  an  equivalent  stress  measure. 
History  effects  are  modeled  with  a  hardening  multiaxial 
formulation  based  on  the  elastic  back  stress  vector  and  an 
equivalent  (scalar)  drag  stress  measure.  Equations  are  derived 
for  an  orthotropic  model  of  incompressible  flow  and  for 
estimating  the  orthotropic  parameters  from  uniaxial  test  data. 

The  uniaxial  model  contains  a  total  of  six  parameters  that 
can  be  determined  from  conventional  experimental  testing  methods 
for  ice.  The  model  is  verified  against  Jacka's  (1984),  Sinha's 
(1978),  and  Brill  and  Camp's  test  data  on  the  creep  of 
polycrystalline  ice.  Predictions  of  the  ratio  of  transient  to 
total  strain  agree  qualitatively  with  Sinha's  equation  if  grain 
size  effects  are  taken  into  account.  The  mechanical  behavior 
under  monotonically  increasing  stress  is  successfully  predicted 
using  Sinha's  (1983)  data  obtained  from  constant  displacement 
rate  tests. 

For  the  multiaxial  model,  experimental  verification  is  made 
difficult  by  the  lack  of  data  for  the  pure  flow  of  freshwater 
polycrystalline  (S2  or  S3)  ice.  The  model  predicts  pressure- 
insensitive  behavior  under  conventional  triaxial  loading 
conditions.  Also,  theoretical  predictions  agree  well  with 
Frederking's  (1977)  data  from  constant  strain  rate  tests  carried 
out  under  plane  strain  conditions  (although  it  should  be  noted 
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that  his  data  is  for  ice  with  distributed  cracks  or  damage 
induced  by  loading,  not  pure  flow.) 
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APPENDIX  A 

Orthotropic  Material  Parameters  (a^^  -  a^ ) 

A. 1  Definition  of  Symbols 

The  derivation  here  considers  only  the  power  law  creep  of 
ice.  The  parameters  a^^  to  ag  and  to  have  been  defined  in 
the  paper.  The  coefficients  a,  b2  to  bg  are  the  constants  for 
the  unaxial  power  law  along  the  y  (reference)-,  x-  and  z- 
directions,  and  along  the  45®  axes  on  the  y-z,  x-y  and  z-x 
planes  respectively.  Thus, 
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(A. 6) 


where  ®yy '••••' ®45 (  2x )  ^yy *^45  (  zx  )  viscous 

strain  rate  and  the  stress  components.  Also,  we  can  set  the 
first  orthotropic  parameter  a^^-l  without  loss  of  generality. 

A. 2  Uniaxial  Tests  in  the  X-  and  Z-  Directions 

To  derive  a2  and  a^,  we  first  obtain  expressions  for 
the  strain  rates  in  the  x-  and  z-directions  using  Eqs.  (50  -  52) 
and  the  definition  for  the  equivalent  stress  (Eq.  45).  Note  that 
the  stress  vectors  are 


—  •  XX 


0  0  0  0  0]^  and  e-tO  0  p,,  0  0  0  )'^ 


for  loading  in  the  X-  and  Z-  directions,  respectively.  Thus; 
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Equations  (A. 7)  and  (A. 8)  are  then  compared  with  Eqs.  (A. 2)  and 

1  /N 

(A. 3),  respectively.  Solving  for  a,  with  <5,-(a/b,)  '  in  the 

^  I  /N  ^ 

first  pair  of  equations  and  with  second  pair 

of  equations,  we  obtain  two  simultaneous  equations  involving  ^2 
and  a^: 

a2  -  -  1  (A. 9) 

a2  -  (a2+a3)e2^*^'^^*^‘^^^  -  1  (A. 10) 

Equations  (69)  and  (70)  are  obtained  by  solving  Eqs.  (A. 9)  and 
(A. 10)  for  a2  and  a^  in  terms  of  and  82* 


A. 3  Uniaxial  Tests  at  45°  on  Y-Z  ,  X-Y  and  Z-X  Planes 

Consider  the  case  of  the  uniaxial  test  at  45®  on  the  y-z 
plane.  The  stress  applied  at  45®  to  the  coordinate  axes  in  the 
plane  is  denoted  by  ‘^45(y2)*  corresponding  strain  is  denoted 

by  c 


45(yz)  • 


By  means  of  a  stress  transformation,  the  plane 


stress  vector  I  1 ’'-I  ,/2  ,  )/2  .  yj  )/2  I 

After  computing  the  equivalent  stress  defined  in  Eq.  (45), 
the  inplane  strains  are  computed  using  Eqs.  (50)>-(52): 


®yy  “  K  a^/6 

®zz  "  ^  ^3^^  ''45(yz)^ 
^yz  “  «5  "45(yz) 


(A. 11) 

(A. 12 ) 

(A. 13) 


where 


K  -  a  (■ 


j(N+l)/2  ^ (1/2)^^"^^(A.14) 


ai+a2 


The  strain  rate  at  45®  to  the  coordinate  axis  can  be  obtained 
by  a  strain  transformation  : 
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ai  +  a3 

Usiyz)-  1/“  ' -  -  2=5)  K  (A.15) 

3 

Comparison  of  Eqs.  (A. 4)  and  (A. 15)  with  a/b^ )  yields  an 
expression  for  a^,  as  given  by  Eq.  (72).  To  obtain  parameters  a^ 
and  ag,  (see  Eqs.  (71)  and  (73))  ,  similar  45*  tests  can  be 
conducted  in  the  x-y  and  z-x  planes  respectively.  For  a 
transversely  isotropic  material, 

constants  a^^  to  a^  can  then  be  simplified  to  the  following: 
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APPENDIX  B 
Frederking's  Tests 


B.l  Type  A  Test 


The  coordinate  axes  are  defined  in  the  text.  The  ice  sheet 
is  subjected  to  normal  stress  in  the  y-direction,  and  its 

in-plane  movement  in  the  z-direction  is  restrained.  Stresses  in 
the  x-direction  are  assumed  to  be  zero.  Thus: 


a  «  0 

XX 


(B.l) 


(B.  2 


The  derivation  below  assumes  that  damage  is  negligible.  Using 
Eqs .  (50)-(52)  and  (B.2),  the  following  expression  is  obtained: 


(B.3) 


After  computing  the  equivalent  stress  (Eq.  (45)),  the  strain 
rate  in  the  y-diiection  is  determined  from  Eqs.  (50)-(52): 


®2®3 
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(B.4) 


where  the  superscript  c  on  implies  that  it  is  confined.  For 

an  unconfined  test  we  have  from  Eq.  (A.l): 


'yy  ■ 


(B.5) 


where  the  superscript  u  on  Cyy  implies  that  it  is  unconfined.  If 
the  strain  rates  are  the  same,  we  can  equate  Eqs.  (B.4)  and 
(B.5)  to  obtain  (with  substitutions  from  Eqs.  ( A. 16 )-( A. 21 )  for 
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transverse  isotropy)  the  expression  for  r  given  by  Eq.  (74) 

z 


The  load  is  applied  in  the  y-direction.  Stresses  in  the 
z-direction  are  assumed  to  be  zero.  Displacements  are  constrained 
in  the  x-direction.  These  imply: 


The  same  procedure  is  followed 
equations  corresponding  to  Eqs. 


(B.6) 

(B.7) 

as  in  the  type  A  test.  The 
(B.3)-(B.4)  are,  respectively: 


(B.8) 


(B.9) 


Comparing  Eqs.  (B.9)  and  (B.5),  and  substituting  from  Eqs. 
(A. 16 )-(A. 21 ) ,  Eq.  (75)  for  follows. 
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APPENDIX  C 
Triaxial  Test 


In  the  triaxial  test  of  a  transversely  isotropic  ice  sheet 

in 

state  is  described  by  the  vector  f  o. 


subjected  to  a  normal  stress  a  in  the  y-direction,  the  stress 

,T 


e  a 

XX  yy  zz 


ace] 
xy  yz  zx ' 


[to  o  to  0  0  0)  ,  where  T  is  the  ratio  of  the  confining  stress 
to  the  axial  stress.  The  equivalent  stress  (Eq.  (45))  is  o^  - 
(1-t)o.  The  strain  rate  in  the  y-direction  is  obtained  from  Eqs. 
(50)-(52)  as  follows: 
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-  a(l-T 


tr  ,N 
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(C.l  ) 


where  the  superscript  'tr'  signifies  loading  under  triaxial 
conditions.  Combining  Eqs.  (B.S)  and  (C.l)  yields  Eq.  (76).  Also 
the  shear  stress  normalized  by  the  unconfined  stress  is 
independent  of  t  as  shown  below: 


(1-T) 
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FIGURE  CAPTIONS 


Figure  1  Dimensionless  Strain  Rate  Plotted  Against 

Dimensionless  Time,  from  the  Data  of  Jacka 
(1984)  . 

Figure  2  Dimensionless  Strain  Rate  Plotted  Against 

Dimensionless  Strain,  from  the  Data  of  Jacka 
(1984  )  . 

Figure  3  Dimensionless  Strain  Plotted  Against 

Dimensionless  Time,  from  the  Data  of  Jacka 
(  1984  )  . 

Figure  4  Dimensionless  Strain  Plotted  Against 

Dimensionless  Time,  from  the  Data  of  Sinha 
(1978)  . 

Figure  5  Strain  Dependence  of  Ratio  of  Transient  Strain 

to  Total  Strain  for  Various  Grain  Sizes; 

0-1. 0  MPa  at  -ICC. 

Figure  6  Stress  and  Strain  History  and  Stress-Strain 

Results  on  Columnar-Grained  Ice  of  Average 
Grain  Diameter  of  4_^5  mm. at  -10®C  for  Nominal 
Strain  Rate  of  5x10”  s”  . 

Comparison  Between  the  Predicted  and  Experimental 
Creep  and  Recovery  of  Snow  Ice.  Experimental  Data 
of  Brill  and  Camp  Reproduced  from  Sinha  (1979). 
The  Zero  Time  of  Curve  (c)  is  Shifted  for 
Clarity. 
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A  RATE-SENSITIVE  CONTINUUM  DAMAGE  MODEL  FOR  ICE 


A  RATE- SENSITIVE  CONTINUUM  OAMAGE  MODEl  FOR  ICE 


Mao  S.  Wu  and  S.  Shyam  Sunder 

Massachusetts  Institute  of  Technology,  Department  of  Civil  Engineering,  Rm 
'1-274,  Cambridge,  KA  C2139  (U.S.A.) 


ABSTRACT 

A  rate-sensitive  damage  model  is  developed  for  describing  the  continuum 
behaviour  of  ice  -  under  variable  loading  conditions.  The  damage  formulation 
employs  a  second  rank  symmetric  tensor  for  describing  the  damage  state  of  the 
material,  a  net  stress  tensor  for  specifying  the  evolution  of  damage,  and  an 
averaged  net  stress  tensor  for  specifying  the  stress-strain  laws.  The  rate  cf 
damage  is  assumed  to  depend  on  an  isotropic  component  associated  with  the 
distortional  stress  and  an  anisotropic  component  associated  with  the  maximum 
principal  deviatoric  stress.  Pressure  melting  is  accounted  for  by  means  of  a 
relation  between  the  creep  resistance  parameter  in  Glen's  power  law  and 
temperature  corrected  for  the  effect  of  hydrostatic  stress.  Integrated  with 
the  generalized  Maxwell  formulation  and  assuming  an  initially  isotropic 
material,  the  model  can  (a)  describe  anisotropic  damage  behaviour  under 
ffiultiaxial  and  non-steady  states  of  stress,  (b)  capture  the  pressure  sensitiv¬ 
ity  cf  ice,  end  (c)  model  pressure  melting  of  ice  under  large  hydrostatic 
stresses.  Verification  cf  the  model  is  achieved  with  several  independent  sets 
of  data,  including  tnose  for  first-year  sea  ice  and  freshwater  ice. 

1 .  INTRODUCTION 

The  theory  of  plasticity  has  been  sucessfully  applied  to  ductile 
materials,  but  is  generally  unsuitable  for  describing  the  behaviour  cf  brittle 
materials.  Tnis  is  because  it  was  originally  developed  as  a  phenomenological 
thecry  for  modelling  nonlinear  beha\'iour  due  to  such  irreversible 
microstructural  changes  as  crystalline  slip  and  twinning,  whereas  the  nonline- 


arity  of  brittle  materials  is  primarily  due  to  the  nucleation  and  growth  of 
microcracks.  It  is  this  inherent  deficiency  in  the  theory  of  plasticity  that 
has  encouraged  the  application  of  damage  mechanics  in  formulating  the 
constitutive  laws  for  a  general  material  in  w'n.ch  several  distinctly  different 
mechanisms,  including  microcracking  and  rate-independent  (or  depende  t)  visco¬ 
us  processes,  may  occur  simultaneously. 

The  concept  of  damage  as  an  internal  variable  which  reflects  the 
accumulation  of  internal  defects  was  pioneered  by  Kachanov  (1956),  vrtio  used 
the  complementary  idea  of  'continuity'  as  a  variable  for  modelling  brittle 
creep  rupture  of  materials.  Since  then,  damage  mechanics  has  been  used  for 
modelling  the  behaviour  of  brittle  rock  (Kachanov,  1982;  Costin,  1963, 
1985a, b),  concrete  (Fonseka  and  Krajcinovic  ,  1981;  Ortiz,  1965;  Wu,  1985), 
and  of  course  metals  (Leckie  and  Hayhurst,  1974;  Murakami  and  Ohno,  1978; 
Chaboche,  1981,  1964;  Z>e%y,  1985). 

Ice  under  compressive  loading  behaves  as  a  continuum  undergoing  damage, 
which  is  particularly  significant  in  the  ductile  to  brittle  transition  region. 
Indeed,  damage  leads  to  strain  softening  under  constant  strainrate  loading, 
and  to  tertiary  creep  under  constant  stress  loading.  Moreover  ,  unloading  a 
damaged  material  often  shows  a  reduction  in  the  elastic  modulus.  These 
observations  nave  stimulated  the  application  of  damage  concepts  in  the 
constitutive  modelling  of  ice  (Karr,  1985;  Szyszkowski  and  Glockner,  1986). 
However,  the  former  model  is  limited  to  one-dimensional  applications,  and  both 
models  assume  that  damage  can  be  described  by  a  scalar  variable,  i.e.,  the 
anisotropic  nature  of  damage  is  ignored.  This  is  in  contradiction  to 
microstructural  observations  of  many  materials  which  show  that  creep  defects 
can  be  orientated  on  planes  perpendicular  to  the  maximum  principal  tensile 
stress.  Such  anisotropy  is  especially  important  under  non-proportional  loading 


when  conditions  for  rupture  or  failure  often  vary  according  to  the  particular 
loading  history  of  the  material. 

Damage  in  ice  is  also  k-.own  to  be  sensitive  to  strainrate,  i.e.,  the 
evolution  of  damage  is  closely  related  to  the  development  of  creep  strain. 
This  is  because  in  the  creep  damage  process  the  nucleation  and  growth  of 
microcavities  can  be  influenced  by  diffusion,  by  dislocation  creep,  and  by 
stress  concentration  at  irregularities  or.  the  grain  boundary  caused  by  grain 
boundary  sliding  (Sinha,  1984).  All  these  processes  are  linked  to  the 
development  of  creep  strain,  which  therefore  affects  damage  growth. 

Ice  is  also  a  pressure  sensitive  material,  particularly  under  high 
strainrate  loading,  which  means  it  has  a  higher  shear  strength  under 
confinement  than  under  uniaxial  compression.  This  is  due  to  the  suppression  of 
iricrocracking  by  the  confining  stress.  For  very  low  strainrate  loading, 
microcracking  is  not  significant  and  little  or  no  pressure  sensitivity  is 
observed.  On  the  other  hand,  it  is  known  that  the  shear  resistance  of  ice 
decreases  under  large  confining  stress,  and  ultimately  vanishes  at  the 
pressure  corresponding  to  phase  change.  This  mode  of  microstructural  change, 
together  with  the  phenomenon  of  pressure  sensitivity  discussed  above,  have 
been  considered  indirectly  through  the  formulation  of  three-dimensional  yield 
envelopes  based  on  experimental  data  of  conventional  triaxial  tests  (Nadreau, 
1986) . 

The  present  study  develops  a  more  elaborate  continuum  damage  theory  for 
the  model  proposed  by  Ting  and  Shyam  Sunder  Cl 985,  1986  ).  It  integrates  an 
anisotropic  damage  theory  and  a  pressure  melting  model  with  the  generalized 
MB:'well  formulation  adopted  in  the  original  model.  Based  on  the  assumption  of 
an  initially  isotropic  material,  the  current  model  can  (a)  descrioe  anisotrop¬ 
ic  damage  behaviour  under  complex  stress  fields,  (b)  capture  the  pressure 


sensitivity  characteristic  of  ice,  and  (c)  model  the  pressure  melting 
behaviour  of  ice  under  large  hydrostatic  stresses.  Verification  of  the  model 
is  achieved  with  several  independent  sets  of  data,  including  those  for  first 
year  sea  ice  and  freshwater  ice. 

2.  CONTINUUM  DAMAGE  MODELLING  OF  ICE 

The  formulation  of  continuum  damage  models  typically  consists  of  three 
critical  steps.  These  are  (a)  the  selection  of  a  damage  variable  based  on  an 
appropriate  averaging  over  a  rs-presentative  volume  ,  (b)  the  establishme;  t  of 
the  damage  evolution  laws,  and  (c)  the  adoption  of  a  set  of  constitutive 
relations.  Intense  debate  exists  as  to  how  a  damage  model  should  be 
constructed.  Krajcinovic  (1983a,  b,  1985a,  b,  c)  firmly  supported  the  use  of  a 
vectorial  damage  variable  to  describe  a  field  of  flat,  planar  microcrachs, 
while  the  traditional  predilection  for  second  rank  tensors  has  led  to  such 
representation  of  the  damage  state  by  Kachanov  (1980),  Murakami  and  Ohno 
(1981),  and  many  others.  The  formulation  of  the  damage  evolution  equations  can 
be  accomplished  by  either  prescribing  an  independent  law  for  each  component  of 
the  damage  tensor  (Chaboche,  1981;  Murakami  and  Ohno,  1981;  Costin,  1985a,  b) 
or  by  using  the  concept  cf  a  damage  surface  (Dragon  and  Mroz ,  1979; 
Krajcinovic  and  Fonseka,  1981  ).  Lack  cf  space  precludes  a  detailed  review  cf 
these  various  approaches;  for  critical  appraisal  see  Krajcinovic  (1985a,  b) 
and  Kachanov  (1985).  Me  shall  discuss  only  the  damage  theory  of  Murakami  and 
Ohno  (1961)  upon  which  the  damage  formulation  is  based. 

Damage  Tensor  :-In  this  theory,  the  microcracks  are  directly  used  to 
define  the  damage  variable  tensor  through  the  reduction  in  the  strength  of 
sections.  For  an  elementary  volume  V  of  the  damaged  material  continuum 
containing  n  microcracks  ,  it  was  shown  that  the  damage  state  can  be  defined 


by  the  following  second  rank  symmetric  tensor  D: 


n 

£  “  -  I  /v  Be  X  Be 

S(V)/3  k*1 


(1) 


where  dSj^  is  the  area  of  the  k-th  microcrack,  Nj^  the  unit  vector  normal  to 
dSjj,  and  S(V)  the  total  area  of  grain  boundaries  in  V.  The  symbol  x  denotes 
tensor  product  of  two  vectors,  the  resulting  quantity  being  a  second  rank 
tensor.  can  be  rewritten  in  terms  of  its  principal  values  Dj  and 
corresponding  principal  directions  Kj  {j  *  1,2,3)  so  that 


where  Dj  represents  the  surface  density  of  the  microcracks  in  the  principal 
plane  perpendicular  to  Kj  (Fig.  la).  The  damage  effects  on  an  arbitrary  plane 
can  be  examined  by  means  of  the  Cauchy  tetrahedron  OABC  (Fig.  lb).  It  shows 
that  due  to  damage  on  the  three  principal  planes,  the  plane  ikBC,  defined  by 
the  undamaged  area  vector  S_v  where  S  is  the  area  of  and  _v  the  unit  normal  to 
the  surface,  can  be  visualized  as  transforming  into  the  fictitious  plane 
A'B'C  described  by  the  net  area  vector  S'^'.  It  can  be  showTi  that  the  two 
area  vectors  are  related  through  the  tensor  (_1_  ; 


S'V'  *  (2  -  D)S_v 


(3) 


Tne  damage  representation  described  above  means  that  (a)  an  arbitrary  damage 


field  can  be  described  in  terms  of  three  mutually  orthogonal  systems  of 
parallel  microcrac)ts,  (b)  the  effect  of  a  multitude  of  damage  fields  is  in 
general  obtained  by  summing  the  different  damage  tensors  (one  for  each  field), 
(c)  the  effect  of  damage  in  a  plane  of  arbitrary  orientation  is  to  transform 
it  into  a  plane  of  diminished  area  and  different  orientation,  and  (d)  the 
tensor  (_1_  ~U)  in  Eq.(3)  can  be  interpreted  as  a  linear  transformation  between 
the  area  vector  S_v  and  the  net  area  vector  S'_v'. 

Damage  Effect  •»nsor  and  Net  Stress  Tensor  :»The  reduction  in  area  of  a 
plane  magnifies  the  stress  tensor.  If  o_  denotes  the  Cauchy  stress  acting  on 
the  tetrahedron  OABC  (Fig. 1b),  then  the  force  vector  acting  on  the  surface  S_v 
is  obtained  using  Cauchy's  formula 

St_-£(S_v)  (4) 

where  z  is  the  corresponding  stress  vector.  The  effective  surface  clement  S'_v' 
is  subjected  to  the  same  force  vector.  Equating  the  two  force  vectors  results 
in 

£'  (S'_v’  )  »  £(££)  (5) 

which  upon  substitution  from  Eq.(3]  gives 

C  (S’ V  )  -  0(1  -  D)-"'  (S'V  )  (6) 


In  Eq.  (5)  above,  c'  is  the  ne 


stress  tensor  and  its  relationship  with  the 


Cauchy  stress  is  evident  from  Eq.  (6),  i.e 


£'  *  -  D)-'>  (7) 

The  net  stress  o'  is  not  necessarily  symmetric?  for  the  purpose  of 
computation  we  use  only  its  symmetric  part 

£'«£(££  +  _*£)  (8) 

where 

£  «(£-£)-■>  (9) 

The  tensor  d  expresses  the  effect  of  damage  on  the  Cauchy  stress  and  is 

termed  the  damage  effect  tensor.  It  can  also  be  thought  of  as  a  quantity 

which  represents  the  effect  of  local  stress  concentration  which  influences  the 
growth  of  damage.  Physically,  the  above  formulation  implies  that  the  effect  of 
the  Cauchy  stress  acting  on  a  volume  element  OABC  of  a  damaged  continuum  is 

identical  to  that  of  the  net  stress  acting  on  the  corresponding  fictitious 

element  O'A'B'C  without  damage.  The  net  stress  is  thus  the  locally  magnified 
stress  which  governs  the  evolution  of  damage,  i.e.,  the  damage  growth  laws  are 
specified  with  respect  to  the  net  stress  tensor.  Eq.  (7)  or  Ec.  (8)  is  the 
three-dimensional  generalization  of  the  net  stress  employed  in  the  classical 
damage  theory  of  Kachanov  and  Rabotnov. 

Net  Stress  Tensor  for  Constitutive  Relations  In  Eq.  (8)  it  is  seen 
that  the  transformation  between  the  Cauchy  stress  tensor  and  the  net  stress 
tensor  is  specified  by  the  second  rank  damage  effect  tensor,  although  a 


trar.sf ormation  between  two  second  rank  tensors  should  be  specified  by  a  fourth 
rank  tensor.  The  result  of  this  transformation  is  that  of  the  two  subscripts 
of  the  Cauchy  stress  tensor,  only  one  is  transformed.  The  implicit  assumption 
is  that  creep  damage  affects  the  stress  only  through  the  change  of  the 
effective  area  on  which  stress  acts,  and  not  through  their  three-dimensional 
arrangement.  While  this  is  generally  true  for  damage  growth,  the  deformation 
characteristics  of  damaged  materials  may  depend  on  both  area  reduction  and  the 
three-dimensional  arrangement  of  the  microdefects.  The  net  stress  tensor 
described  above  therefore  strictly  cannot  be  used  to  specify  the  constitutive 
relations . 

Murakami  and  co-workers  (1982,  1983,  1986  )  proposed  the  use  of  a  fou-th 
rank  net  stress  tensor  constructed  from  a  rather  complex  fourtn  rank  damage 
effect  tensor  .  However,  by  measuring  the  uniaxial  tensile  stress  of  brass 
specimens  perforated  with  holes  at  different  angles,  they  found  that  if  the 
cavity  density  was  below  a  certain  limit,  the  uniaxial  tensile  stress  was 
essentially  the  same  regardless  of  the  orientations  of  the  perforations,  and 
that  this  limit  usually  exceeds  that  corresponding  to  rupture  of  the  material. 
This  implies  that  the  effect  of  material  damage  on  creep  deformation  may  be 
isotropic  when  the  cavity  density  is  less  than  a  critical  limit,  and  -.is  may 
be  valid  ^ust  oefore  failure.  In  view  of  these  experimental  observations,  it 
was  proposed  that  the  net  stress  tensor  for  constitutive  rela cions  (averaged 
net  stress  tensor)  may  be  defined  through  a  scalar  invariant  of  the  damage 
effect  tensor  as  follows 

-  (n  +  (1  -  ri)/3  -r  £  )£  0  <  n  <  1  (10) 


wnere  r, 


is  a  material  constant. 


Evolution  Ecuations  of  Damaae: 


must  be  stated  at  the  outset  that  it 


is  the  net  stress  that  is  used  in  the  damage  laws.  To  avoid  cumbersome 
terminology,  however,  we  shall  drop  the  word  'net*  in  the  description  of  the 
stresses.  For  the  proposed  model,  it  is  postulated  that  the  rate  of  damage 
growth  is  governed  by  the  local  state  of  stress  ,  the  effect  of  the  local 
stress  concentration  at  the  microcavities  expressed  through  the  damage  effect 
tensor  the  effective  total  strainrate  Cg  and  the  temperature  T,  i.e., 

£  “  »  1'  ^e'  ")  ( ^ 

To  obtain  Ec.  ("I)  in  explicit  form,  we  can  maXe  use  of  tensor  function 
theory  as  well  as  experimental  observations.  For  many  materials  iricrocracks 
may  form  in  planes  perpendicular  to  the  principal  tensile  stress  ,  and  may 
also  develop  isotropically.  Under  a  non-hydrostatic  compressive  state  of 
stress,  local  tensile  fields  can  arise  from  material  property  mismatch  between 
grains  or  from  contact  stresses  between  grains  with  irregular  grain  boundaries 
(Costin,  1963).  Tnese  local  tensile  stresses  will  result  in  the  nucleation  and 
growth  cf  microcracks.  Since  the  normals  to  these  mhcrocracks  are  usually 
orientated  in  tne  directions  cf  the  principal  tensile  deviatoric  stresses,  it 
can  be  assumed  that  the  local  tensile  stresses  are  proportional  in  magnitude 
to  the  principal  tensile  deviatoric  stresses  {henceforth  principal  deviatoric 
stresses)  emd  alsc  act  in  the  same  directions  .  Thus,  damage  may  be  assumed  to 
develop  in  the  directions  ^  and  (i  »  1,  2,  3)  of  the  positive  (tensile) 
principal  values  of  £'  and  its  deviator  _S ' ,  as  well  as  isotropically  in  all 
three  principal  directions.  For  this  model,  we  further  postulate  that  damage 
develops  in  the  direction  of  the  maximum  principal  deviatoric  stress  since  it 


is  predominantly  in  this  direction  that  microcracXs  evolve  under  axial 
compression.  This  is  referred  to  as  the  principal  stress  damage  law  (Murakami 
and  Ohno,  1978). 

For  many  brittle  materials,  the  volumetric  or  bulk  stress  is  known  to 
suppress  the  formation  and  growth  of  microcracks.  This  phenomenon  is  also  true 
for  ice,  especially  at  high  rates  of  loading  (Mellor,  1963),  and  can  be 
modelled  by  mitigating  the  damage  due  to  the  maximum  principal  deviatoric 
stress.  The  amount  of  damage  that  is  suppressed  depends  on  the  magnitude  of 
the  volumetric  stress.  For  very  large  volumetric  stress,  the  damage  due  to  the 
maximum  principal  deviatoric  stress  r^ay  be  completely  suppressed.  Seen  from  a 
different  angle,  the  above  formulation  is  equivalent  to  a  decouplir ;  cf  the 
applied  stresses  into  deviatoric  and  volumetric  components  with  different 
weighting  factors. 

It  is  conceivable  that  material  damage  can  also  be  attributed  to  the 
cistortional  stress,  or  the  effective  stress.  This  isotropic  damage  is  related 
to  the  second  invariant  cf  the  deviatoric  stress  tensor,  and  the  damage 
evolution  is  termed  the  J2  damage  law  (Murakami  and  Ohno,  1978)  .  The 
volumetric  stress  suppresses  the  microcracks  and  leads  to  increase  in  material 
strength  and  hence  the  cistortional  stress.  The  large  cistortional  stress  will 
induce  damage  and  consequently  lim.its  the  gain  in  strength  due  to  the  increase 
in  volumetric  stress.  Note  that  the  J2  damage  law  may  become  dominant  under 
moderate  and  large  confining  stresses  when  anisotropic  damage  is  significantly 
suppressed  . 

The  above  reasoning  implies  that  the  damage  evolution  laws  consist  cf 
two  parts  -an  anisotropic  part  due  to  the  maximum  principal  deviatoric  stress, 
appropriately  reduced  according  to  the  magnitude  of  the  hydrostatic  stress, 
and  an  isotropic  part  due  to  the  effective  stress. 


To  incorporate  the  effect  of  strainrate  on  damage,  it  is  proposed  that  a 
power  law  holds  between  the  damage  growth  rate  and  the  effective  total 
strainrate,  i.e,, 

D  a  (12) 

The  effect  of  the  local  stress  concentration  on  damage  evolution  can  be  ta)ten 
into  account  by  incorporating  invariants  of  the  damage  effect  tensor  in  the 
damage  evolution  equations.  For  tne  isotropic  part,  the  invariant  (1/3 
tr(t2)1/2j  which  is  related  to  the  root  mean  square  of  the  diagonal  terms  of 
the  damage  effect  tensor,  is  adopted  .  For  the  anisotropic  part,  the  invariant 
tr(c  X  )  is  used,  where  denotes  the  direction  cosine  vector  of  the 
maximum  principal  deviatoric  stress.  This  expression  simply  evaluates  the 
magnitude  of  the  damage  effect  titsso^  in  the  direction  of  the  maximum 
principal  deviatoric  stress.  Thus,  for  the  case  of  plane  stress,  »  (cos 
s'*,  sin  eM,  where  S'*  is  the  angle  between  the  maximum  principal  deviatoric 
stress  and  the  reference  x-axis,  and  tr(_^  x  )  is  the  familiar 
expression  cos^  6^  +  ^22  sin^  6^. 

The  principal  damage  axes  and  general  do  not  remain 
fixed  with  respect  to  the  reference  axes  x^ ,  X2  and  X3,  especially  for  the 
case  of  a  non-steady  state  of  stress  in  which  the  maximum  principal  deviatoric 
stress  may  rotate  in  an  arbitrary  fashion  %irith  time.  The  relative  amount  of 
damage  in  the  xi ,  X2,  X3,  XTX2f  X2''*3  *3"Xl  directions  at  any  instant  of 
time  is  determined  by  the  tensor  product  x  ^1,  For  plane  stress  the 
relative  amount  of  damage  is  cos^  S'*  in  the  x^  direction,  sin^  6^  in  the  X2 
direction,  and  sin  e"*  cos  e"*  in  the  xtX2  direction.  I 


f  there  are  two  equal 


principal  stresses,  as  in  the  case  of  conventional  triaxial  loading,  we  have 
to  sum  the  contributions  from  each  principal  stress. 

The  proposed  damage  evolution  equations  are  based  on  the  concepts 
discussed  above.  They  can  be  written  in  the  following  form  : 

£  “  iA  +  il  (13) 

where  A  and  I  denote  'anisotropic'  and  'isotropic'  respectively,  and 

^  •=  £3!  ^  £3''  <051  '  +  Bl,  '/3>N  X  (14) 

il  *  Ce=1  (1/3  tr£2)C2/2  (^ ( 3j2« ) 1/2)N  £  (15) 

In  Ecs.(14)  and  (15),  ci ,  03 »  Of  6  and  y  are  material  constants  .  N  is  the 
power  law  index  for  ice;  here  we  have  assumed  that  the  same  power  governs 
the  deformation  characteristics  as  well  as  the  damage  behaviour.  S',  li'/3 
and  (3^2*)'''^^  denote  the  maximum  principal  deviatoric  stress,  the  volumetric 
stress  and  the  effective  stress,  respectively.  The  cash  sign  indicates  net 
stress,  i.e.. 


si' 

*  max ( '),  i  ■  1,  2 1 

3 

(16) 

-1’ 

-  Ii(o') 

(17) 

^2' 

-  J2(ff’) 

(18) 

and  the 

quantity  inside  the 

angle 

brackets,  if  negative,  is  assigned  a  value 

of  zero.  Other  symbols  are  as  defined  previously. 


Constitutive  Relations: 


-The  constitutive  laws  for  the  proposed  model 


are  based  on  the  generalized  Maxwell  formulation  adopted  in  the  original  model 
of  Ting  and  Shyam  Sunder  (’985).  These  laws  are  reproduced  below  for  the 
isotropic  case;  they  are  merely  a  special  case  of  the  original  model  with  the 
orthot.  jpic  parameters  a<f  to  ag  equal  to  1 .  It  must  be  emphasized  that 
'stress'  in  this  section  actually  refers  to  the  'averaged  net  stress'  defined 
by  Eq.  (10). 

The  total  strainrate  tensor  e_  can  be  decomposed  into  the  elastic 

a  • 

strainrate  tensor  _C  £  and  the  creep  strainrate  tensor  The  latter  can  be 
further  decoupled  into  the  delayed  elastic  tensor  and  the  secondary  creep 
strainrate  tenser.  Damage  is  incorporated  through  the  averaged  net  stress 
tensor,  and  affects  all  the  strain  components  identified  above,  thus  resulting 
in  tertiary  creep.  The  governing  equations,  in  engineering  notation  of 
stresses  and  strains,  are  : 


•  •  • 

t  =  C  C-  +  E, 


£cr  *  ''“a 


(19) 

(20) 


wnere 


-I  i  'S'  *  ‘X  > 


(21  ) 


In  Eq.  (20),  is  the  conventional  deviatoric  stress  tensor,  in  Ec.  (21) 

e 

is  the  effective  st.ess  for  isotropic  materials.  The  two  terms  in  Eq.  (21)  are 
associated  with  the  delayed  elastic  component  and  the  secondary  creep  strain 
component,  the  latter  being  the  multiaxial  form  of  Glen's  power  law.  M  and  A 


are  creep  resistance  parameters.  is  the  modulus  of  elasticity  given  by 


by 


•  E{^/z  exp(A/E  -  1) 


(22) 


where  E  is  Young's  modulus  and  r  a  constant. 

Pending  experimental  verification,  we  assume  in  this  model  that  the 
averaged  net  stress  is  given  by  Eq.  (10)  with  n  *  0,  that  is, 

-  (1/3  tr^)£  (23) 

This  implies  that  the  net  stress  affecting  the  global  stress-strain  behaviour 
is  governed  only  by  the  average  damage  effect  on  stress. 

3.  PRESSURE  MELTING  MODEL 

Ice  is  )tnown  to  change  into  liquid  form  at  a  certain  level  of 
hydrostatic  pressure,  'biaxial  testing  of  ice  was  carried  out  only  relatively 
recently  (Jones,  1978,  1982  ;  Hausler,  1981;  Richter-Menge  et  al.  ,  1985? 
Nadreau,  1986),  but  analyses  have  already  beer,  done  on  the  available  data  to 
determine  a  three-dimensional  yield  envelope  for  tee.  Unfortunately,  most  of 
these  analyses  do  not  taXe  into  account  the  decreasing  resistance  of  ice  as 
the  hydrostatic  pressure  increases. 

Nadreau  (1986)  proposed  a  new  yield  envelope  which  was  calibrated  with 
his  own  test  data  and  also  results  performed  on  freshwater  ice  from  Jones. 
This  three-dimensional  envelope  is  teardrop-shaped  and  for  isotropic  ice  is 
symmetrical  about  the  hydrostatic  axis.  The  site  of  this  envelope  increases  as 
the  strainrate  increases.  Indeed,  according  to  Mellor  (1983),  for  very  small 


strainrates  the  deviatoric  yield  stress  decreases  steadily  as  the  hydrostatic 
stress  increases  until  it  vanishes  at  the  pressure  corresponding  to  phase 
change.  The  yield  envelope  is  a  cor.e  symmetrical  about  the  hydrostatic  axis. 
For  high  strainrates,  the  deviatoric  yield  stress  increases  initially,  but 
eventually  decreases  as  the  hydrostatic  stress  becomes  sufficiently  large.  The 
yield  envelope  is  thus  similar  to  the  teardrop  surface  proposed  by  Nadreau. 

The  current  model  will  predict  increasing  deviatoric  yield  stress  with 
increasing  hydrostatic  stress.  It  will  also  predict  a  peak  with  subsequent 
reduction  in  the  deviatoric  yield  stress.  This  occurs  when  the  hydrostatic 
stress  IS  large  enough  so  that  any  suppression  of  damage  is  offset  by  the 
increase  in  damage  due  to  the  large  distortional  stress.  At  large  hydrostatic 
stress,  however,  the  phenomenon  of  pressure  melting  becomes  important.  In 
order  to  incorporate  this  effect,  it  is  proposed  that  the  creep  resistance 
parameter  A  in  Eq.  (2i)  decreases  with  increasing  hydrostatic  stress.  The 
pressure  melting  model  will  be  briefly  discussed  in  the  following  paragraphs. 

The  value  cf  the  creep  resistance  parameter  A  at  any  temperature  T  and 
hydrostatic  pressure  p  is  obtaimd  through  a  temperature  correction,  that  is, 
is  equal  to  the  value  at  an  equivalent  temperature  T'  and  zero  hydrostatic 
pressure.  Mathematically,  this  is  expressed  as 


A(T  ,  p)  -  A(T'  ,  0) 

where  the  equivalent  temperature  is  defined  as 


(24) 


(To  -  Tir.) 


(25) 


In  Eq.  (25),  Tq  is  the  melting  temperature  at  zero  hydrostatic  pressure,  and 


Tm  is  the  melting  temperature  at  some  hydrostatic  pressure  p.  The  value  of 
can  be  ofctained  from  the  phase  diagram  shown  in  Fig. 2  (from  Nadreau,  1986)  , 
or  from  the  following  relation  proposed  for  ice  1h,  valid  between  To=  273  K 
for  pure  ice  (or  ''71  K  for  sea  ice)  and  some  reference  temperature  T.j  (e.g. 
251  K),  and  between  0  MPa  and  some  reference  hydrostatic  pressure  p^  (e.g.  200 
MPc )  ; 


-mtP)  -  +  (To  -  TO  (1  -  (p/Pi)^)’/’  (26) 

where  i  is  a  constant  which  describes  the  curvature  of  the  temperature  versus 
pressure  curve.  Typrcally,  ^  varies  between  1  and  1.2;  for  the  former  value, 
the  curve  is  linear.  An  average  value  of  1.05  is  used  in  this  study.  Also, 
the  trend  of  experimental  data  (Richter-Menge  et  al.,  1985)  suggests  that  an 
appropriate  value  for  p^  to  be  50  MPa  for  sea  ice. 

Having  obtained  the  value  of  the  equivalent  temperature  using  Eqs.(26) 
and  (25),  we  can  determine  the  value  of  P  according  to  the  following  relations 

A  =  Ac  ex?[  2/(N?"')  1  T'  C  263  K  (27) 

A/A(2e3  ,  C)  *  :  (To  -  T’)/(To  -  263)  T'  >  263  K  (28) 

where  A(263  ,  0)  is  the  creep  resistance  at  T  ■  263  K  (and  p  «  0  MPa)  and  ;  is 
a  constant  determined  by  the  requirement  of  continuity  in  the  slope  between 
the  two  curves  at  T  *  263  K  (see  Fig.  3),  C  is  equal  to  0.377  for  pure  ice  and 
equal  to  0.301  for  sea  ice  for  this  formulation.  Q,  N  and  R  are  the  activation 
energy,  the  power  law  index  for  ice  and  the  universal  gas  constant  respective¬ 
ly.  Aq  is  a  temperature  independent  constant. 


4.  EXPERIMENTAL  VERIFICATION 


The  model  is  verified  against  several  independent  sets  of  experimental 

data  obtained  from  constant  strainrate  tests.  The  following  values  f^r  the 

model  parameters  provide  an  average  representation  of  the  data: 

E  -  9500  MPa 

r  -  0.65 

Aq  -  9.2  MPa  sV^’ 

N  «  3 

M  -  1411 

Q  m  65000  J  mol-'' 

The  universal  gas  constant  is  equal  to  8.314  C  mol”^  K-\  The  damage  model 
introduces  five  material  dependent  parameters  a,  6,  1,  and  C2  to  model 

anisctropic/isctropic  damage,  pressure  sensitivity,  strainrate  sensitivity  and 
the  effect  of  local  stress  concentration  on  damage  .  More  precisely,  a,  6  and 
T  are  associated  with  the  maximum  principal  deviatoric  stress,  the  volumetric 
stress  and  the  effective  stress,  respectively  .  Ci  is  a  constant  which 
determines  the  rate  cf  damage  due  to  the  strainrate  effect  while  C2  is 
associated  with  the  contribution  to  damage  growth  due  tc  the  damage  effect  or 
the  effect  of  local  stress  concentration. 

Typically,  a  lies  between  1  and  1C,  6  between  0  and  20,  and  y  between  0 
and  1.  Tneir  units  are  given  by  MPa”^  s^-l”’'^/^'.  Note  the  rather  wide  range  cf 
values  o  can  assume:  a  large  a  vrili  result  in  the  lowering  of  the  maximum 

stress  versus  strainrate  curve  in  the  ductile  to  brittle  transition  region 
(case  1 )  ,  while  a  smaller  o  will  cause  departure  from  the  power  law  model 
prediction,  although  not  necessarily  result  in  a  lowering  of  the  above- 
mentioned  curve  (case  2)  .a  also  determines  the  magnitude  of  the  peak  stress. 
The  ratio  a/B  is  usually  greater  than  about  2.  y  is  usually  small  compared  to 
a  ,  reflecting  the  fact  that  damage  evolution  is  generally  governed  by  the 
principal  stress  law  rather  than  the  J2  law.  However,  the  J2  law  can  be 


important  under  large  hydrostatic  stresses  when  the  damage  due  to  the  maximum 
principal  stress  is  significantly  suppressed.  This  suggests  that  the  material 
deforms  in  a  pseudo-ductile  manner  with  large  scale  deformation  being  taken  up 
by  many  homogeneously  distributed  cracks,  as  is  the  case  for  brittle  solids 
with  microcracks  or  pores  under  large  compressive  stress  states  (Ashby 

and  H?llam  ,  1986).  c^  typically  lies  between  0.2  and  2.5,  while  C2  can 

assume  any  value  between  0  and  2.  The  smaller  the  value  of  c^ ,  the  greater 

will  be  the  damage  rate.  If  C2  is  large,  very  severe  post-peak  strain 

softening  under  constant  strainrate  loading  can  be  expected.  Table  1  gives 
the  values  of  the  five  parameters  which  pro%'ide  an  average  representation  of 
available  experimental  data  corresponding  to  zero  salinity  and  T  «  -10®C  for 

the  two  cases  mentioned  above. 

Fig.  4  shows  the  maximum  stress  or  strength  observed  from  constant 
strainrate  tests  on  columnar  pure  and  sea  ice  versus  strainrate  (case  1).  The 
experimental  data  has  been  normalized  for  the  temperature  of  -10®C  and  zero 
salinity.  Tne  sclid  line  in  the  figure  represents  the  prediction  cf  the 
current  model  using  the  first  set  of  damage  parameters  in  Table  1  ;  it  is  seen 
that  the  model  captures  the  overall  trend  cf  the  data  very  well,  particularly 
'the  severe  strain-softening  behaviour  in  the  ductile  to  brittle  transition 
region.  For  strainrates  greater  than  10”^  s”"'  the  continuum  model  is  invalid 
and  a  horizontal  line  representing  sudden  failure  by  brittle  fracture  is  drawn 
at  the  stress  level  of  5  MPa.  Tne  dashed  line  indicates  the  prediction  of  the 
power  law  model  which  does  not  consider  in  its  formulation  the  concept  of 
damage.  Ihe  effect  of  temperature  on  the  strength  versus  strainrate  behaviour 
is  shown  in  Fig.  5,  which  shows  that  w’ith  decreasing  temperature  the  maximum 
strength  increases  and  also  occurs  at  lower  strainrates.  At  the  strainrates 
of  1C"^  and  10"^  the  strength  at  -40*C  is  about  twice  and  quadruple  that 


at  -10®C,  respectively.  The  same  plot  also  suggests  that  in  the  region 
dominated  by  brittle  behaviour,  the  strength  becomes  insensitive  to  temperatu¬ 
re  .  These  predictions  agree  well  with  Mellor's  observations  (see  Fig.  35  cf 
Mellor,  1983).  The  evolution  of  damage  parallel  to  the  compressive  axis  with 
respect  to  strains  and  strainrates  is  illustrated  in  Fig.,  6.  It  is  observed 
that  damage  is  insignificant  at  very  small  strains  and  strainrates,  and 
eventually  levels  off.  Tne  latter  observation  implies  t-iat  the  evolution  of 
axial  damage  under  compression  is  relatively  stable  compared  to  the  tensile 
fracture  of  a  brittle  solid  (Ashby  and  Hallam,  1986;  Nemat-Nasser  and  Horii, 
1982) . 

Some  researchers  (e.g.  Mellor,  1963)  have  suggested  that  the  maximum 
stress  versus  strainrate  behaviour  deviates  from  the  power  law  predictions 
but  does  not  exhibit  a  peak  with  subsequent  reduction,  as  illustrated  in  Fig, 
7.  This  corresponds  to  case  (2)  in  the  present  study.  The  model  predictions 
using  the  second  set  of  damage  parameters  in  Table  1  and  15  *  2.5  are 
superimposed  on  the  experimental  data.  The  motivation  for  us-nc  a  smaller  N 
in  this  single  instance  is  to  match  the  slope  of  the  curve  at  very  small 
strainrates,  but  this  also  requires  the  use  of  a  larger  A,  equal  to  650000 
.MTa  £  '  It  is  seen  tl-iat  the  model  can  capture  the  characteristics  displayed 
by  Gcld's  data  for  various  ice  types  corresponding  to  the  temperature  cf 
-10*C  and  density  of  0.9  Mg/m^. 

The  proposed  model  has  also  been  verified  against  the  triaxial  test 
results  of  Richter-Menge  et  al.  (1985),  Jones  (1982)  and  Nadreau  (1986). 

In  the  tests  performed  by  the  former,  the  columnar  first-year  sea  ice 
samples  were  horizontally  cored  in  three  directions  and  at  different  depths  of 
the  ice  block.  The  direction  cf  testing  is  specified  by  the  angle  between  the 
compressive  axis  and  the  c-axis  of  the  ice  (c  :  c).  The  confining  pressure  is 


ramped  in  fixed  proportion  to  the  applied  stress  in  the  uniaxial  direction. 
The  four  ratios  of  confining  pressure  to  axial  pressure  used  in  their  tests 
are  7  —  Of  0.25f  0.5  and  0.75. 

The  model  predictions  of  the  stress-strain  behaviour  for  at  c  *  0°  and 
for  the  two  strainrates  of  10“3  and  10“S  g-l  aj-g  shown  in  Fig.  6.  The  same 
figure  also  plots  the  experimental  data,  suitably  averaged  through  the 
thickness  of  the  ice  block.  Note  the  suppression  of  damage  reflected  through 
the  decrease  in  strain-softening  in  the  stress-strain  behaviour  as  the 
confining  pressure  increases.  The  corresponding  normalized  shear  strength 
versus  confining  pressure  cvrves  are  plotted  in  Fig.  9.  It  is  observed  that 
the  pressure  sensitivity  of  sea  ice  is  moderate,  with  only  about  25%  increase 
in  shear  strength.  As  7  equals  0.75,  pressure  melting  has  become  dominant.  It 
can  be  concluded  that  the  model  gives  reasonbly  good  predictions  cf  the 
characteristics  of  triaxial  behaviour,  that  is,  pressure  sensitivity  and 
pressure  melting,  at  least  for  the  strainrate  of  1C"^  s"'*.  For  the  lower 
strainrate  cf  10"^  s“''.  Fig.  9  shows  widely  scattered  data  and  this  explains 
why  the  corresponding  stress-strain  curves  in  Fig.  8  do  not  fit  closely  with 
the  data  for  this  strainrate. 

The  same  general  conilusions  can  be  drawn  with  respect  to  the  other  two 
cases  cf  c  :  c  =  45°  and  90°  ,  see  Figs.  10  and  11.  From  the  shear  strength 
versus  confining  pressure  plots  it  is  observed  that  pressure  sensitivity  is 
most  significant  when  the  loading  axis  is  aligned  with  the  c-axis,  and  least 
significant  when  they  are  at  an  angle  of  90°  (only  15%  increase  in  shear 
strength ) . 

The  proposed  model  assumes  that  the  material  is  initially  isotropic, 
while  most  of  the  sea  ice  samples  that  were  tested  are  materially  anisotropic. 
To  overcome  this  problem,  verification  of  the  model  has  been  carried  out  by 


assuming  that  the  creep  resistance  A  in  different  directions  is  different,  and 
also  that  the  damage  parameters,  being  material  dependent  constants,  also  vary 
with  the  direction  of  loading.  In  other  words,  if  the  material  were 
isotropic,  A  and  the  damage  parameters  would  be  invariant  with  respect  to  the 
direction  of  loading.  For  the  model  A(0*^)/A(  45°)  ■  1.7,  and  A{0°)/A(90°)  = 
1.2.  It  was  also  found  that  for  sea  ice  6  and  y  are  much  larger  than  the 
corresponding  values  for  non-saline  ice.  Tne  first  three  rows  of  Table  2  give 
the  values  of  these  parameters. 

The  effect  of  confining  pressure  on  damage  parallel  to  the  compressive 
axis  is  illustrated  in  Fig.  12.  The  suppression  of  damage  due  to  the  confining 
pressure  is  very  evident.  Tne  decomposition  of  damage  into  its  anisotropic  and 
isotropic  components  is  shown  in  Fig.  13.  In  Fig.  13a  it  is  seen  that  the 
anisotropic  damage  due  to  tne  maximum  principal  deviatoric  stress  is  complete¬ 
ly  suppressed  for  t  >  C.S,  while  Fig.  13b  shows  that  the  accompanying  increase 
in  the  distortional  stress  due  to  the  increase  in  confining  stress  actually 
aggravates  damage.  However,  at  i:  =  C.75,  pressure  melting  becomes  dominant  and 
the  distortional  stress  that  can  be  sustained  by  the  ice  reduces  significantly 
,  thus  resulting  in  smaller  damage,  as  seen  in  Fig.  13b.  The  crossing-over  cf 
the  curves  at  small  strains  in  Fig. 13b  is  attributed  to  the  fact  that 
initially  the  distortional  stress  and  hence  the  isotropic  damage  reduces  with 
increase  in  confining  stress j  this  is  overridden  when  the  enhanced  suppression 
of  anisotropic  damage  by  the  larger  confining  stress  leads  to  increase  in  the 
distortional  stress.  Finally,  it  should  be  mentioned  that  since  the  maximum 
principal  deviatoric  stress  is  always  normal  to  the  compressive  axis  and  has 
no  component  along  this  axis,  the  damage  in  planes  perpendicular  to  the 
compressive  axis  is  solely  due  to  the  distortional  stress  and  is  represented 


by  Fig.  13b 


Jones  (1982)  performed  a  series  of  conventional  triaxial  tests  on 
freshwater  ice  at  the  temperature  of  -12“C  and  for  strainrates  varying  between 
1.4  X  10~6  £-1  and  1.4  x  10"^  s“\  Fig.  14  reproduces  his  data  and  also  plots 
the  model  predictions  using  the  damage  parameters  listed  in  Table  2.  The  model 
gives  good  predictions,  particularly  for  strainrates  at  and  under  1.4  x  10“2 
s"'';  for  the  very  high  strainrates  (5.4  x  10"3  and  1.4  x  10“^  ,  the  model 
slightly  overpredicts  the  shear  strength.  Note  that  the  shear  strength  can 
increase  by  as  much  as  100%  under  confinement  for  the  highest  strainrate 
used  in  his  tests. 

To  formulate  a  three-dimensional  yield  envelope  for  ice  which  takes  into 
account  pressure  sensitivity  and  pressure  melting,  Nadreau  has  recently 
conducted  a  series  of  conventional  triaxial  tests  on  saline  ice  and  iceberg 
ice.  Both  types  of  ice  are  isotropic.  The  two  strainrates  used  in  his  tests 
are  10~^  and  1C”®  s"’’.  The  experimental  results  for  saline  ice  at  the 
temperatures  of  263  K  and  268  K,  and  for  iceberg  ice  at  263  K,  are  shown  in 
Figs.  15  and  16  respectively  .  Superimposed  on  these  data  are  the  predictions 
of  the  model,  which  can  be  seen  to  capture  the  overall  behaviour  very  well. 
For  the  smaller  strainrate,  the  model  predicts  little  pressure  sensitivity. 
The  iceberg  ice  also  displays  a  greater  degree  of  pressure  sensitivity  than 
the  saltne  ice,  with  approximately  100%  and  65%  increase  in  shear  strength 
respectively.  The  values  of  A  and  the  damage  parameters  are  shown  in  Table  2. 

5.  CONCLUSIONS 

Tne  model  presented  in  this  paper  for  describing  the  continuum  damage 
behaviour  of  ice  under  variable  loading  conditions  is  based  on  the  generalized 
Maxwell  differential  formulation  ,  an  anisotropic  damage  theory  and  a  pressure 
melting  model.  Specifically,  it  is  able  to  (a)  describe  anisotropic  damage 


behaviour  under  mulriaxial  and  non-steady  states  of  stress,  (b)  capture  the 
pressure  sensitivity  of  ice,  and  (c)  model  the  pressure  melting  behaviour  of 
ice.  Vrrif ication  of  the  model  is  achieved  with  several  independent  sets  of 
data,  including  thos  for  first-year  sea  ice  and  freshwater  ice.  The  following 
conclusions  can  be  drawn  from  the  work  reported  in  this  paper: 

1.  The  damage  model  is  described  by  5  parameters,  which  are  used  in  the 
damage  evolution  equation;  to  model  the  dependence  of  damage  growth  rate  on 
the  principal  stress,  the  volumetric  stress,  the  effective  stress,  as  well  as 
strainrate  and  the  effect  of  local  damage. 

2.  Tne  damage  model  assumes  that  (a)  the  damage  state  can  be  described 
by  a  second  rank  symmetric  tensor,  (b)  the  net  stress  tensor,  which  is  the 
locally  magnified  stress  due  to  area  reduction,  can  be  used  to  specify  the 
damage  growth  laws,  and  (c)  an  averaged  net  stress  tensor  is  valid  for 
specifying  the  constitutive  laws. 

3.  Material  damage  is  considered  to  consist  of  an  isotropic  component 
and  an  anisotropic  component.  The  former  is  associated  with  the  J2  damage  law 
while  the  latter  with  the  principal  stress  damage  law.  Damage  due  to  the 
distortional  stress  may  be  dominant  under  moderate  or  large  hydrostatic 
stresses  when  anisotropic  damage  is  suppressed. 

4.  The  pressure  melting  model  is  based  on  the  assumption  that  the  creep 
resistance  parameter  A  reduces  with  increase  in  hydrostatic  pressure.  This 
is  mathematically  achieved  through  a  temperature  correction  using  the  phase 
diagram  for  ice  or  a  mathematical  equation.  The  relation  between  A  and  the 
corrected  temperature  is  then  obtained  via  the  Arrhenius  Law  for  T  <  263  K  or 
a  newly  formulated  relation  for  T  >  263  K. 

5.  The  model  can  capture  well  the  dominant  characteristics  cf  the 
experimental  data  by  Wang,  Mellor,  Richter-Menge  et  al.,  Jones,  Nadreau  and 


t 

others  .  In  the  uniaxial  tests,  the  behaviour  of  decreasing  strength  or 
asymptotically  increasing  strenqth  with  increasing  strainrate  in  the  brittle 
to  ductile  transition  region  can  be  modelled.  At  strainrates  greater  than  10”^ 
s"'',  the  presence  of  macrocracks  precludes  a  solely  continuum  description  of 
ice  behaviour.  Damage  is  insignificant  at  small  strains  and  strainrates. 
Also,  strength  becomes  insensitive  to  temperature  at  the  brittle  end  of  the 
strength  versus  strainrate  curve.  For  the  triaxial  tests,  the  effect  of 
microcrack  suppression  by  the  volumetric  stress  can  be  captured  by  the  model  , 
as  verified  against  the  stress-strain  curves  of  Richter-Menge  et  al.  Pressure 
sensitivity  iz  moderate  for  sea  ice,  but  is  much  more  significant  for 
freshwater  ice.  Under  large  hydrostatic  stresses,  pressure  melting  becomes 
important. 

Additional  research  is  needed  to  resolve  several  issues;  including  (a) 
the  experimental  verification  of  the  damage  tensor,  the  damage  effect  tensor, 
and  the  net  stress  tensors  for  damage  evolution  and  constitutive  relations, 
(b)  the  development  of  the  model  assuming  initial  material  anisotropy  as  well 
as  initial  distribution  of  pcres/damage,  (c)  a  more  in-depth  exploration  of 
the  micromechanics  of  damage  such  as  the  possibility  of  complete  or  partial 
crack  closure  due  to  rotation  of  the  stress  field  with  possible  effect  on 
damage  accumulation,  the  interaction  of  microcracks  leading  to  instability  , 
and  the  influence  of  the  shape  of  the  microdefects  (planar/spheroidal/ 
ellipsoidal)  on  material  damage  behaviour,  (d)  the  application  of  the  model  to 
stress  controlled  tests  and  arbitrary  loading  histories,  and  (e)  further 
investigation  into  the  relationship  between  creep  resistance  and  the  hydrosta¬ 
tic  stress. 
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Table  1  :  Average  Values  of  A  and  Damage  Parameters  (for  Zero  Salinity  and 
tne  Temperature  of  -10®C) 
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BM 
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Case  2 

1.45 

0.04 

0.0001 

0.1 

650000 

Table  2  :  Specific  Values  of  A  and  Damage  Parameters 

a  6  y  c-i  C2  A 

Ricnter-Menge  et.  al 
(Columnar  Sea  Ice) 
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1  .05 
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0.31 
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a  :  C  *  45® 

1  .59 

0.36 

0.07 

0.31 
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C  ;  C  =  90® 

1.10 

C.50 

0.06 

0.31 

0.1 

127500 

Jones 

(Freshwater  Ice) 
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(Saline  Ice) 
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